
import cartopy
import cartopy.crs as ccrs
import collections
import datetime
import gsw
import gudinfo
import IPython.display
import netcdf_tools
import math
import matplotlib as mpl
import matplotlib.colors
import matplotlib.image
import matplotlib.patches
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import matplotlib.transforms
import numpy as np
import os
import os.path
import pandas as pd
import PIL
import PIL.Image
import read_mitgcm
import scipy
import scipy.interpolate
from scipy.io import loadmat
import vstats


plt.close("all")


'''latlon2xy convenience function to convert latlon to chosen projection,
and figure object'''
def latlon2xy(lon,lat,proj):
    xyz = proj.transform_points(ccrs.PlateCarree(),lon,lat).squeeze()
    x,y,_ = np.split(xyz,3,axis=len(xyz.shape)-1)
    x[np.isinf(x)]=np.nan
    y[np.isinf(y)]=np.nan
    return x.squeeze(),y.squeeze()


# STATIONS
Station=collections.namedtuple(
    'Station','lon lat marker  markeredgecolor markerfacecolor markersize')
station1=Station(-60.3919999999999,68.3096,'o','none','red',5)
station2=Station(-61.463,68.0358,'o','none','red',5)
station3=Station(-62.355,67.8668,'o','none','red',5)
station4=Station(-63.5639999999999,67.54,'o','none','red',5)
station5=Station(-63.783,67.482,'*','black','yellow',10) # Qikiqtarjuaq SIC
station6=Station(-64.638,67.239,'o','none','red',5)
stations=[station1,station2,station3,station4,station5,station6]


plotbathymetry=True  

with plt.style.context('mplstyles/maps.mplstyle'):
    fig=plt.figure(figsize=(10, 16))
    
    # --- CCGS AMUNDSEN 2018, LEG 2B
    proj_am=ccrs.LambertConformal(central_longitude=-63,
                                  central_latitude=67.5,
                                  standard_parallels=(67,68))
    ax_am=fig.add_subplot(223,projection = proj_am)
    ax_am.coastlines(resolution='10m')
    ax_am.add_feature(cartopy.feature.LAND)
    if not plotbathymetry:
        ax_am.add_feature(cartopy.feature.OCEAN)
    MapExtent = [-66,-60,66,69]
    ax_am.set_extent(MapExtent, ccrs.PlateCarree())
    gl = ax_am.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                         linewidth=1, color='gray', alpha=0.5,
                         linestyle='--')
    gl.xlocator = mticker.FixedLocator(np.arange(-68,-58))
    gl.ylocator = mticker.FixedLocator(np.arange(65,70))

    # stations
    for station in stations:
        xstation,ystation = latlon2xy(np.array([station.lon]),
                                      np.array([station.lat]),
                                      proj_am)
        marker=station.marker
        markerfacecolor=station.markerfacecolor
        markersize=station.markersize
        markeredgecolor=station.markeredgecolor
        ax_am.plot(xstation,ystation,marker,markeredgecolor=markeredgecolor,
                   markerfacecolor=markerfacecolor,markersize=markersize)
        
    # longitude labels
    y0,_ = ax_am.get_ylim()
    for lon in [-66,-64,-62,-60]:
        x0,_ = proj_am.transform_point(lon,MapExtent[2]-1,
                                       src_crs = ccrs.PlateCarree())
        plt.text(x0,y0-2e4,'{:2d}$^\circ$W'.format(-lon),
                 horizontalalignment = 'center',verticalalignment='top')

    # latitude labels
    x0,_ = ax_am.get_xlim()
    someLons = np.arange(-80,-50)
    for lat in np.arange(66,69+1):
        # interpolate latitude circle to map boundary
        xyzArray = proj_am.transform_points(ccrs.PlateCarree(),
                                            someLons,lat*np.ones_like(someLons))
        x = xyzArray[:,0]
        y = xyzArray[:,1]
        y0 = np.interp(x0,x,y)
        plt.text(x0-2e4,y0,'{:2d}$^\circ$N'.format(lat),
                 horizontalalignment = 'right',verticalalignment='center')
    
    if plotbathymetry:
        # color bar for bathymetry
        cmap = plt.cm.get_cmap('Greys')
        cmap = cmap.from_list('Custom cmap',
                              [cmap(i) for i in range(150,50,-1)], 100)
        cmap.set_over('white')
        bounds = [0,500,1000,2000]
        norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

    # IBCAO bathymetry
    if plotbathymetry:
        top = loadmat('data/DataS1_observations_IBCAO_1min_bathy.mat')
        lon_ibcao,lat_ibcao=np.meshgrid(top['x'],top['y'])
        z_ibcao=-top['IBCAO_1min']
        hb_am=ax_am.contourf(*latlon2xy(lon_ibcao,lat_ibcao,proj_am),
                             z_ibcao,levels=[0,500,1000,2000],
                             cmap=cmap,norm=norm,extend='max')
        axcb_am= fig.colorbar(hb_am,ax=ax_am,orientation='horizontal')
        axcb_am.ax.set_title('Water depth (m)')
        axcb_am.ax.set_position([0.70,-0.10, 0.25,0.20])
    
    # labels for stations
    for i,station in enumerate(stations):
        xstation,ystation = latlon2xy(np.array([station.lon]),
                                      np.array([station.lat]),
                                      proj_am)
        ax_am.text(xstation+3e3,ystation+1e4,i+1,
                   horizontalalignment = 'center',
                   verticalalignment='center')
    # text
    b_lon,b_lat = -66.5, 69.1,
    ax_am.text(b_lon,b_lat,'(b)',
               horizontalalignment = 'left',fontsize=30,
               transform=ccrs.PlateCarree())
    
    # --- BAFFIN BAY
    proj_bb=ccrs.LambertConformal(central_longitude=-60,
                                  central_latitude=70,
                                  standard_parallels=(65,75))
    ax_bb=fig.add_subplot(221,projection=proj_bb)
    ax_bb.coastlines(resolution='50m')
    ax_bb.add_feature(cartopy.feature.LAND)
    if not plotbathymetry:
        ax_bb.add_feature(cartopy.feature.OCEAN)
    MapExtent = [-76,-48,62,79]
    ax_bb.set_extent(MapExtent, ccrs.PlateCarree())
    gl = ax_bb.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                         linewidth=1, color='gray', alpha=0.5,
                         linestyle='--')
    gl.xlocator = mticker.FixedLocator([-80,-75,-70,-65,-60,
                                        -55,-50,-45,-40,-35])
    gl.ylocator = mticker.FixedLocator(np.arange(58,84,3))

    # stations
    for station in stations:
        xstation,ystation = latlon2xy(np.array([station.lon]),
                                      np.array([station.lat]),
                                      proj_bb)
        marker=station.marker
        markerfacecolor=station.markerfacecolor
        markersize=station.markersize
        markeredgecolor=station.markeredgecolor
        ax_bb.plot(xstation,ystation,marker,markeredgecolor=markeredgecolor,
                   markerfacecolor=markerfacecolor,markersize=markersize)

    # longitude labels
    y0,_ = ax_bb.get_ylim()
    for lon in [-70,-60,-50]:
        x0,_ = proj_bb.transform_point(lon,MapExtent[2]-1,
                                       src_crs = ccrs.PlateCarree())
        plt.text(x0,y0-2e4,'{:2d}$^\circ$W'.format(-lon),
                 horizontalalignment = 'center',verticalalignment='top',
                fontsize=20)

    # latitude labels
    x0,_ = ax_bb.get_xlim()
    someLons = np.arange(-120,-40,1)
    for lat in np.arange(64,78,3):
        # interpolate latitude circle to map boundary
        xyzArray = proj_bb.transform_points(ccrs.PlateCarree(),
                                            someLons,lat*np.ones_like(someLons))
        x = xyzArray[:,0]
        y = xyzArray[:,1]
        y0 = np.interp(x0,x,y)
        plt.text(x0-2e4,y0,'{:2d}$^\circ$N'.format(lat),
                 horizontalalignment = 'right',verticalalignment='center')
    
    # IBCAO bathymetry
    if plotbathymetry:
        hb_bb=ax_bb.contourf(*latlon2xy(lon_ibcao,lat_ibcao,proj_bb),
                             z_ibcao,levels=[0,500,1000,2000],
                             cmap=cmap,norm=norm,extend='max')

    # text
    a_lon,a_lat = -98, 77.5
    ax_bb.text(a_lon,a_lat,'(a)',
               horizontalalignment = 'left',fontsize=30,
               transform=ccrs.PlateCarree())
    greenland_lon,greenland_lat = -38, 77
    ax_bb.text(greenland_lon,greenland_lat,'Greenland',
               horizontalalignment = 'right',
               transform=ccrs.PlateCarree())
    baffinbay_lon,baffinbay_lat = -60, 72
    ax_bb.text(baffinbay_lon,baffinbay_lat,'Baffin\nBay',
               horizontalalignment = 'right',
               transform=ccrs.PlateCarree())

    # rectangle
    extentStudyarea = ax_am.get_extent()
    vertices = np.array([
        [extentStudyarea[0], extentStudyarea[2]], # left, bottom
        [extentStudyarea[0], extentStudyarea[3]], # left, top
        [extentStudyarea[1], extentStudyarea[3]], # right, top
        [extentStudyarea[1], extentStudyarea[2]], # right, bottom
        [extentStudyarea[0], extentStudyarea[2]]
    ])
    xVert,yVert = vertices.transpose()
    xBox = np.interp(np.arange(0,4,.02),np.arange(0,4.1,1),xVert)
    yBox = np.interp(np.arange(0,4,.02),np.arange(0,4.1,1),yVert)
    xyzBoxLatlon = ccrs.PlateCarree().transform_points(ax_am.projection,
                                                       xBox,yBox)
    xBoxLatlon = xyzBoxLatlon[:,0]
    yBoxLatlon = xyzBoxLatlon[:,1]
    xyzBoxTransformed = proj_bb.transform_points(ccrs.PlateCarree(),
                                                 xBoxLatlon,yBoxLatlon)
    xBoxTransformed = xyzBoxTransformed[:,0]
    yBoxTransformed = xyzBoxTransformed[:,1]
    ax_bb.plot(xBoxTransformed,yBoxTransformed,'r')
    ax_bb.relim()
    
    # --- WORLD MAP
    proj_wm = ccrs.NearsidePerspective(central_longitude=-60,
                                       central_latitude=60)
    ax_wm = fig.add_subplot(222,projection = proj_wm)
    ax_wm.set_global()
    ax_wm.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='black')
    ax_wm.add_feature(cartopy.feature.OCEAN)
    gl = ax_wm.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                         linewidth=1, color='gray', alpha=0.5,
                         linestyle='--')
    extentStudyarea = ax_bb.get_extent()
    vertices = np.array([
        [extentStudyarea[0], extentStudyarea[2]], # left, bottom
        [extentStudyarea[0], extentStudyarea[3]], # left, top
        [extentStudyarea[1], extentStudyarea[3]], # right, top
        [extentStudyarea[1], extentStudyarea[2]], # right, bottom
        [extentStudyarea[0], extentStudyarea[2]]
    ])
    xVert,yVert = vertices.transpose()
    xBox = np.interp(np.arange(0,4,.02),np.arange(0,4.1,1),xVert)
    yBox = np.interp(np.arange(0,4,.02),np.arange(0,4.1,1),yVert)
    xyzBoxLatlon = ccrs.PlateCarree().transform_points(ax_bb.projection,
                                                       xBox,yBox)
    xBoxLatlon = xyzBoxLatlon[:,0]
    yBoxLatlon = xyzBoxLatlon[:,1]
    xyzBoxTransformed = proj_wm.transform_points(ccrs.PlateCarree(),
                                                 xBoxLatlon,yBoxLatlon)
    xBoxTransformed = xyzBoxTransformed[:,0]
    yBoxTransformed = xyzBoxTransformed[:,1]
    ax_wm.plot(xBoxTransformed,yBoxTransformed,'r')
    
    ax_wm.set_position( [0.60,0.60 ,0.38,0.38])
    ax_bb.set_position( [0.10,0.53 ,0.60,0.42])
    ax_am.set_position( [0.10,0.05 ,0.60,0.42])

    # --- SAVE
    if plotbathymetry:
        plt.savefig('figures_progress/maps.png')
    else:
        plt.savefig('figures_progress/maps_nobathy.png')


years=range(1,11)
nbyears=len(years)
iyeartempo=nbyears
first_year=np.arange(0,365)
last_year=np.arange(365*(nbyears-1)+0,365*(nbyears-1)+365)
array1d_iT1y_iT=last_year
first_year366=np.arange(0,366) # first year for heatmaps


gridfile='data/DataS8_output_mitgcm/exp0/grid.t001.nc';
drF=netcdf_tools.read_netcdf(gridfile,'drF')


RC=-netcdf_tools.read_netcdf(gridfile,'RC')


RF=-netcdf_tools.read_netcdf(gridfile,'RF')


RF_above151=np.array(RF[RF<151])


obsice=np.empty(365)
obsice[:]=np.NaN
obsice[0:199]=1
obsice[199]=0


obs_breakup=next(i for i,v in enumerate(obsice) if v < 0.50)


def load_thickness():
    fname='data/DataS2_observations_Qikiqtarjuaq/66407.csv'
    df = (
        pd.read_csv(
            fname,
            dtype={
                'mission': 'string',
                'date': str,
                'date_time': str,
                'latitude': np.float32,
                'longitude': np.float32,
                'sample_type': 'category',
                'sample_thickness_cm_average': np.float32,
                'sample_thickness_cm_sd': np.float32,
                'pi': 'category'
            },
            parse_dates=['date','date_time']
        )
    )
    df['doy']=df['date'].apply(lambda x:x.timetuple().tm_yday)
    df['sample_thickness_m_average']=df['sample_thickness_cm_average'] \
    .apply(lambda x:x/100)
    return df


thickness_df=load_thickness()


snow_thickness_df=thickness_df[(thickness_df.mission=='ice_camp_2016') \
                               & (thickness_df.sample_type=='snow')]


ice_thickness_df=thickness_df[(thickness_df.mission=='ice_camp_2016') \
                              & (thickness_df.sample_type=='ice')]


def load_par():
    fname='data/DataS2_observations_Qikiqtarjuaq/66399.csv'
    df = (
        pd.read_csv(
            fname,
            dtype={
                'mission': 'string',
                'station': 'category',
                'date': str,
                'snow_thickness': 'category',
                'depth_m': np.float32,
                'pi': 'category',
                'averaged_par_d_fit_muein_m_2_s_1': np.float32,
                'averaged_par_d_fit_percent_percent': np.float32,
                'averaged_par_d_fit_daily_ein_m_2_day_1': np.float32,
                'averaged_par_d_noon1hloc_ein_m_2_1h_1': np.float32,
                'averaged_par_d_p1h_ein_m_2_1h_1': np.float32,
                'averaged_par_d_p3h_ein_m_2_3h_1': np.float32,
                'averaged_par_d_p24h_ein_m_2_day_1': np.float32,
                'averaged_par_d_p48h_ein_m_2_48h_1': np.float32
            },
            parse_dates=['date']
        )
        .rename(
            columns=dict(
                depth_m='depth',
                averaged_par_d_fit_daily_ein_m_2_day_1='par',
            )
        )
    )
    df['doy']=df['date'].apply(lambda x:x.timetuple().tm_yday)
    return df


par_df=load_par()


par_df=par_df[par_df.mission=='ice_camp_2016']


par_df['par'][par_df['par']<0]

Series([], Name: par, dtype: float32)


g=par_df.groupby('doy')['depth'].value_counts()
g.where(g>1).dropna()

Series([], Name: depth, dtype: float64)


parwater_df=(par_df.dropna(subset=['par'])).copy()
for doy in pd.unique(parwater_df['doy']):
    ice_thickness=parwater_df.loc[(parwater_df.doy==doy),('depth')].min() # m
    parwater_df.loc[parwater_df.doy==doy,'depth']\
    =parwater_df.loc[parwater_df.doy==doy,'depth']-ice_thickness


array1d_iT_obsPAR3m=np.empty(365)
array1d_iT_obsPAR3m[:]=np.NaN
for doy in pd.unique(parwater_df['doy']):
    y_interp=scipy.interpolate.interp1d(
        parwater_df[parwater_df.doy==doy]['depth'],
        parwater_df[parwater_df.doy==doy]['par'])
    array1d_iT_obsPAR3m[doy-1]=y_interp(3)


array1d_iT_obsisolume=np.empty(365)
array1d_iT_obsisolume[:]=np.NaN
c_isolume=0.415 # mol photons m^-2 d^-1
for doy in pd.unique(parwater_df['doy']):
    parwater_onedoy_df=parwater_df[parwater_df.doy==doy]
    if parwater_onedoy_df.iloc[0,:]['par'] > c_isolume:
        idx_isolume=vstats.find_idx_nearest \
        (array=parwater_onedoy_df['par'],
         value=c_isolume)
        isolume=parwater_onedoy_df['depth'].iloc[idx_isolume]
        array1d_iT_obsisolume[doy-1]=isolume


def load_nutrients():
    fname='data/DataS2_observations_Qikiqtarjuaq/66412.csv'
    df = (
        pd.read_csv(
            fname,
            dtype={
                'mission': 'string',
                'date': str,
                'date_time': str,
                'latitude': np.float32,
                'longitude': np.float32,
                'sample_type': 'category',
                'sample_source': 'category',
                'snow_thickness': 'category',
                'code_operation': 'category',
                'station': 'category',
                'station_type': 'category',
                'cast': np.float32,
                'bottle': np.float32,
                'bottom_depth_m': np.float32,
                'target_depth_m': 'string',
                'depth_m': np.float32,
                'filter_type': 'category',
                'quality_flag': 'category',
                'dilution_factor': np.float32,
                'no3_um_l': np.float32,
                'no2_um_l': np.float32,
                'no2_and_no3': np.float32,
                'po4_um_l': np.float32,
                'sioh4_um_l': np.float32,
                'n_ice_samples_ice_camp_2016': np.float32,
                'method': 'category',
                'pi': 'category'
            },
            parse_dates=['date','date_time']
        )
        .rename(
            columns=dict(
                depth_m='depth',
                no3_um_l='no3',
                po4_um_l='po4',
                sioh4_um_l='sioh4'
            )
        )
    )
    df['doy']=df['date'].apply(lambda x:x.timetuple().tm_yday)
    return df


nutrients_df=load_nutrients()


nutrients_df=nutrients_df[(nutrients_df.mission=='ice_camp_2016') \
              & (nutrients_df.sample_type=='water') \
             & (nutrients_df.filter_type=='gff')]


g=nutrients_df.groupby('doy')['depth'].value_counts()
g.where(g>1).dropna()

Series([], Name: depth, dtype: float64)


array1d_iT_obsnitrate3m=np.empty(365)
array1d_iT_obsnitrate3m[:]=np.NaN
array1d_iT_obssilica3m=np.empty(365)
array1d_iT_obssilica3m[:]=np.NaN
for doy in pd.unique(nutrients_df['doy']):
    y_interp=scipy.interpolate.interp1d(
        nutrients_df[nutrients_df.doy==doy]['depth'],
        nutrients_df[nutrients_df.doy==doy]['no3'])
    array1d_iT_obsnitrate3m[doy-1]=y_interp(3)
for doy in pd.unique(nutrients_df['doy']):
    y_interp=scipy.interpolate.interp1d(
        nutrients_df[nutrients_df.doy==doy]['depth'],
        nutrients_df[nutrients_df.doy==doy]['sioh4'])
    array1d_iT_obssilica3m[doy-1]=y_interp(3)


fname='data/DataS5_observations_davis_strait.csv'
fet_df=pd.read_csv(
    fname,
    dtype={
        'doy':np.int64,
        'depth':np.float32,
        'Fe_D_CONC_BOTTLE_nmol_kg':np.float32
    }
)
fet_df['depth']=-fet_df['depth']
# fet_df


# nmol kg^-1 ===> mmol m^-3
fet_df['Fe_D_CONC_BOTTLE_nmol_kg']=fet_df['Fe_D_CONC_BOTTLE_nmol_kg']/1000
fet_df.rename({'Fe_D_CONC_BOTTLE_nmol_kg':'Fe_D_CONC_BOTTLE_mmol_m3'},
             axis='columns',
             inplace=True)
# fet_df


def load_pigment():
    fname='data/DataS2_observations_Qikiqtarjuaq/66417.csv'
    df = (
        pd.read_csv(
            fname,
            dtype={
                'mission': 'string',
                'date': str,
                'date_time': str,
                'latitude': np.float32,
                'longitude': np.float32,
                'sample_id': 'string',
                'sample_code': 'string',
                'sample_type': 'category',
                'sample_source': 'category',
                'snow_thickness': 'category',
                'leg': 'category',
                'code_operation': 'category',
                'station': 'category',
                'station_type': 'category',
                'cast': np.float32,
                'bottle': np.float32,
                'target_depth_m': 'string',
                'depth_m': np.float32,
                'pigment': 'category',
                'conc_mg_m2': np.float32,
                'conc_ug_l': np.float32,
                'conc_mg_m3': np.float32,
                'filtered_vol_l': np.float32,
                'vol_after_melt_l': np.float32,
                'ice_added_vol_of_water': np.float32,
                'ice_total_vol_water': np.float32,
                'ice_calculated_vol_ice': np.float32,
                'dilution_factor': np.float32,
                'quality_control_peaks': 'category',
                'nb_core': np.float32,
                'nb_core_9cm_dia': np.float32,
                'core_area_9cm_dia': np.float32,
                'nb_core_14cm_dia': np.float32,
                'core_area_14cm_dia': np.float32,
                'area_m2': np.float32,
                'description_of_sample': 'category',
                'n': np.float32,
                'observations': 'category',
                'method': 'category',
                'pi': 'category'
            },
            parse_dates=['date','date_time']
        )
        .rename(
            columns=dict(
                depth_m='depth',
                conc_ug_l='chlfluo',
                conc_mg_m3='chlHPLC'
            )
        )
    )
    df['doy']=df['date'].apply(lambda x:x.timetuple().tm_yday)
    return df


chlHPLC_df=load_pigment()


chlHPLC_df=chlHPLC_df[(chlHPLC_df.mission=='ice_camp_2016') \
              & (chlHPLC_df.sample_type=='water') \
             & (chlHPLC_df.pigment=='Total Chlorophyll a') \
             & (chlHPLC_df.method=='HPLC')]


chlHPLC_df['chlHPLC'][ chlHPLC_df['chlHPLC']<0 ]

Series([], Name: chlHPLC, dtype: float32)


g=chlHPLC_df.groupby('doy')['depth'].value_counts()
g.where(g>1).dropna()

doy  depth
146  1.5      2.0
     10.0     2.0
Name: depth, dtype: float64


chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC','doy']][(chlHPLC_df.doy==146) | (chlHPLC_df.doy==148)]
chlHPLC_oneday_df


chlHPLC_df.loc[[20364,20366],'doy']=148


chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC','doy']][(chlHPLC_df.doy==146) | (chlHPLC_df.doy==148)]
chlHPLC_oneday_df


g=chlHPLC_df.groupby('doy')['doy'].value_counts()
g.where(g<=2).dropna()

doy  doy
120  120    2.0
Name: doy, dtype: float64


chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC','doy']][(chlHPLC_df.doy==120)]
chlHPLC_oneday_df


chlHPLC_df=chlHPLC_df.drop(chlHPLC_df[(chlHPLC_df.doy==120)].index)


#doy: 123
depths_case1 =np.array([      1.5, 5,   10,   20,    40])
weights_case1=np.array([      3.25,4.25, 7.5, 15,    10])
#doys: 125,127
depths_case2 =np.array([0.5, 1.5, 5,          20,    40])
weights_case2=np.array([1,   2.25,9.25,       17.5,  10])
#doys:130,132,134,137,139,144,146,148,151,153,155,158,160,162,165,167,169,172,176
depths_case3 =np.array([0.5,1.5,  5,   10,   20,    40])
weights_case3=np.array([1,  2.25, 4.25, 7.5, 15,    10])
#doy: 141
depths_case4 =np.array([0.5,1.5,  5,   10,           40])
weights_case4=np.array([1,  2.25, 4.25,17.5,         15])
#doy: 174
depths_case5 =np.array([0.5,1.5,        10,   20,    40])
weights_case5=np.array([1,  4.75,        9.25,15,    10])
#doys: 179,181,183,186,188,190
depths_case6 =np.array([0.5,1.5,        10,   20,30,    60])
weights_case6=np.array([1,  4.75,        9.25,10,15,     0])
#doys: 193,195,200
depths_case7 =np.array([     1.5,        10,   20,30,    60,  75])
weights_case7=np.array([     5.75,        9.25,10,15,     0,   0])
#doys: 197,204
depths_case8 =np.array([     1.5,        10,   20,30,  45, 60])
weights_case8=np.array([     5.75,        9.25,10,12.5, 2.5,0])


nT=365
array1d_iT_obsvint40mchlHPLC=np.full(nT,np.NaN)
for iT in np.arange(0,nT):
    obsvintchlHPLC=np.NaN
    chlHPLC=[]
    doy=iT+1
    chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC']][(chlHPLC_df.doy==doy)]
    depths=chlHPLC_oneday_df['depth']
    if set(depths_case1).issubset(set(depths)):
        for depth in depths_case1:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case1)
    elif set(depths_case2).issubset(set(depths)):
        for depth in depths_case2:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case2)
    elif set(depths_case3).issubset(set(depths)):
        for depth in depths_case3:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case3)
    elif set(depths_case4).issubset(set(depths)):
        for depth in depths_case4:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case4)
    elif set(depths_case5).issubset(set(depths)):
        for depth in depths_case5:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case5)
    elif set(depths_case6).issubset(set(depths)):
        for depth in depths_case6:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case6)
    elif set(depths_case7).issubset(set(depths)):
        for depth in depths_case7:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case7)
    elif set(depths_case8).issubset(set(depths)):
        for depth in depths_case8:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case8)
    array1d_iT_obsvint40mchlHPLC[iT]=obsvintchlHPLC


#doy: 123
depths_case1 =np.array([      1.5, 5,   10,   20,    40])
weights_case1=np.array([      3.25,4.25, 7.5, 15,    70])
#doys: 125,127
depths_case2 =np.array([0.5, 1.5, 5,          20,    40])
weights_case2=np.array([1,   2.25,9.25,       17.5,  70])
#doys:130,132,134,137,139,144,146,148,151,153,155,158,160,162,165,167,169,172,176
depths_case3 =np.array([0.5,1.5,  5,   10,   20,    40])
weights_case3=np.array([1,  2.25, 4.25, 7.5, 15,    70])
#doy: 141
depths_case4 =np.array([0.5,1.5,  5,   10,           40])
weights_case4=np.array([1,  2.25, 4.25,17.5,         75])
#doy: 174
depths_case5 =np.array([0.5,1.5,        10,   20,    40])
weights_case5=np.array([1,  4.75,        9.25,15,    70])
#doys: 179,181,183,186,188,190
depths_case6 =np.array([0.5,1.5,        10,   20,30,    60])
weights_case6=np.array([1,  4.75,        9.25,10,20,    55])
#doys: 193,195,200
depths_case7 =np.array([     1.5,        10,   20,30,    60,  75])
weights_case7=np.array([     5.75,        9.25,10,20,    22.5,32.5])
#doys: 197,204
depths_case8 =np.array([     1.5,        10,   20,30,  45,60])
weights_case8=np.array([     5.75,        9.25,10,12.5,15,47.5])


nT=365
array1d_iT_obsvintchlHPLC=np.full(nT,np.NaN)
for iT in np.arange(0,nT):
    obsvintchlHPLC=np.NaN
    chlHPLC=[]
    doy=iT+1
    chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC']][(chlHPLC_df.doy==doy)]
    depths=chlHPLC_oneday_df['depth']
    if set(depths_case1).issubset(set(depths)):
        for depth in depths_case1:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case1)
    elif set(depths_case2).issubset(set(depths)):
        for depth in depths_case2:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case2)
    elif set(depths_case3).issubset(set(depths)):
        for depth in depths_case3:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case3)
    elif set(depths_case4).issubset(set(depths)):
        for depth in depths_case4:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case4)
    elif set(depths_case5).issubset(set(depths)):
        for depth in depths_case5:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case5)
    elif set(depths_case6).issubset(set(depths)):
        for depth in depths_case6:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case6)
    elif set(depths_case7).issubset(set(depths)):
        for depth in depths_case7:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case7)
    elif set(depths_case8).issubset(set(depths)):
        for depth in depths_case8:
            chlHPLCtempo \
            =chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
            assert chlHPLCtempo.size==1
            chlHPLCtempo=chlHPLCtempo.iat[0]
            chlHPLC.append(chlHPLCtempo)
        obsvintchlHPLC=np.dot(chlHPLC,weights_case8)
    array1d_iT_obsvintchlHPLC[iT]=obsvintchlHPLC


obsphotofile='data/DataS7_literature.csv'
obsphoto_df=pd.read_csv(obsphotofile)


obsPCmax_df=obsphoto_df[obsphoto_df.parameter=='PCmax']


obsalphachl_df=obsphoto_df[obsphoto_df.parameter=='alphastar']


obsEk_df=obsphoto_df[obsphoto_df.parameter=='Ek']


ice=np.fromfile('data/DataS4_output_nemo_lim3/siarea.nemo.2016.365.32bits.bin',
                dtype='>f')


sim_freezeup=next(i for i,v in reversed(list(enumerate(ice))) if v < 0.50)


ice_file='data/DataS4_output_nemo_lim3/Ice_d.nc'


array1d_iT_snow_volumefull=\
netcdf_tools.read_netcdf(ice_file,'snow volume')
array1d_iT_snow_volumefull=array1d_iT_snow_volumefull[:,1,1]


array1d_iT_snow_volume=array1d_iT_snow_volumefull[365*3:365*4]


array1d_iT_ice_volumefull=\
netcdf_tools.read_netcdf(ice_file,'ice volume')
array1d_iT_ice_volumefull=array1d_iT_ice_volumefull[:,1,1]


array1d_iT_ice_volume=array1d_iT_ice_volumefull[365*3:365*4]


indir='data/DataS8_output_mitgcm/exp0'


info=gudinfo.RunInfo(indir)


esd=info.f['esd'].to_numpy()


stod=86400 # seconds to day
pcmax=info.f['pcmax'].to_numpy() # s^-1
pcmax=pcmax*stod                 # d^-1


ksatno3=info.f['ksatno3'].to_numpy()


stod=86400 # seconds to day
gmax=info.f['grazemax'].to_numpy() # s^-1
gmax=gmax*stod                     # d^-1


carfile='data/DataS8_output_mitgcm/exp0/car.0000000000.t001.nc'


stod=86400 # 3600*24 s per d
array2d_idepth_iT_PARfullumolEonm2s=netcdf_tools.read_netcdf(carfile, 'PAR').squeeze().transpose()
if np.all(array2d_idepth_iT_PARfullumolEonm2s==0):
    #define GUD_PARUICE
    array2d_idepth_iT_PAR_icefullumolEonm2s=netcdf_tools.read_netcdf(carfile, 'PAR_ice').squeeze().transpose()
    array2d_idepth_iT_PAR_icefull=array2d_idepth_iT_PAR_icefullumolEonm2s*1e-6*stod
    array2d_idepth_iT_PAR_owfullumolEonm2s=netcdf_tools.read_netcdf(carfile, 'PAR_ow').squeeze().transpose()
    array2d_idepth_iT_PAR_owfull=array2d_idepth_iT_PAR_owfullumolEonm2s*1e-6*stod
    (ndepths,ndays)=array2d_idepth_iT_PAR_icefull.shape
    array2d_idepth_iT_icefull=np.tile(ice,(ndepths,ndays//365))
    array2d_idepth_iT_meanPARfull=array2d_idepth_iT_PAR_owfull*(1-array2d_idepth_iT_icefull)\
    +array2d_idepth_iT_PAR_icefull*array2d_idepth_iT_icefull
else:
    #undef GUD_PARUICE
    array2d_idepth_iT_meanPARfull=array2d_idepth_iT_PARfullumolEonm2s*1e-6*stod
array2d_idepth_iT_meanPAR=array2d_idepth_iT_meanPARfull[:,array1d_iT1y_iT]
(ndepths,ndays)=array2d_idepth_iT_meanPAR.shape
array1d_iT_modPAR3m=np.zeros(ndays)
x=3
xp=[RC[2],RC[3]]
assert np.all(np.diff(xp)>0), "xp is in decreasing order: %r %r" % (xp[0],xp[1])
for iT in first_year:
    fp =[array2d_idepth_iT_meanPAR[3,iT],array2d_idepth_iT_meanPAR[2,iT]]
    PAR3m=np.interp(x,xp,fp)
    array1d_iT_modPAR3m[iT]=PAR3m


array2d_idepth_iT_meanPARumolEonm2s=array2d_idepth_iT_meanPAR/1e-6/stod


c_isolume=0.415 # mol photons m^-2 d^-1
(ndepths,ndays)=array2d_idepth_iT_meanPAR.shape
array1d_iT_modisolume=np.empty(ndays)
array1d_iT_modisolume[:]=np.NaN
for iT in first_year:
    array1d_idepth_meanPAR=array2d_idepth_iT_meanPAR[:,iT]
    if array1d_idepth_meanPAR[0]>c_isolume:
        idx_isolume=vstats.find_idx_nearest(array=array1d_idepth_meanPAR,
                                           value=c_isolume)
        isolume=RC[idx_isolume]
        array1d_iT_modisolume[iT]=isolume


array2d_idepth_iT_modno3full=netcdf_tools.read_netcdf(carfile, 'TRAC04').squeeze().transpose()
array2d_idepth_iT_modno3=array2d_idepth_iT_modno3full[:,array1d_iT1y_iT]
array2d_idepth_iT_modsioh4full=netcdf_tools.read_netcdf(carfile, 'TRAC06').squeeze().transpose()
array2d_idepth_iT_modsioh4=array2d_idepth_iT_modsioh4full[:,array1d_iT1y_iT]
array2d_idepth_iT_modpo4full=netcdf_tools.read_netcdf(carfile, 'TRAC05').squeeze().transpose()
array2d_idepth_iT_modpo4=array2d_idepth_iT_modpo4full[:,array1d_iT1y_iT]


(ndepths,ndays)=array2d_idepth_iT_modno3.shape
array1d_iT_modnitrate3m=np.zeros(ndays)
x=3
xp=[RC[2],RC[3]]
assert np.all(np.diff(xp)>0), "xp is in decreasing order: %r %r" % (xp[0],xp[1])
for iT in first_year:
    fp =[array2d_idepth_iT_modno3[3,iT],array2d_idepth_iT_modno3[2,iT]]
    nitrate3m=np.interp(x,xp,fp)
    array1d_iT_modnitrate3m[iT]=nitrate3m
(ndepths,ndays)=array2d_idepth_iT_modsioh4.shape
array1d_iT_silicicacid3m=np.zeros(ndays)
x=3
xp=[RC[2],RC[3]]
assert np.all(np.diff(xp)>0), "xp is in decreasing order: %r %r" % (xp[0],xp[1])
for iT in first_year:
    fp =[array2d_idepth_iT_modsioh4[3,iT],array2d_idepth_iT_modsioh4[2,iT]]
    silica3m=np.interp(x,xp,fp)
    array1d_iT_silicicacid3m[iT]=silica3m


array2d_idepth_iT_modfetfull=netcdf_tools.read_netcdf(carfile, 'TRAC07').squeeze().transpose()
array2d_idepth_iT_modfet=array2d_idepth_iT_modfetfull[:,array1d_iT1y_iT]


indir_exp0='data/DataS8_output_mitgcm/exp0'


chlfile=os.path.join(indir,'chl.0000000000.t001.nc')
array2d_idepth_iT_modchlfull=\
read_mitgcm.get_array2d_idepth_iT_chlfull(chlfile)
array2d_idepth_iT_modchl=array2d_idepth_iT_modchlfull[:,array1d_iT1y_iT]


# get vertically integrated Chl a for one year}
def get_array1d_iT_vintchl(indir,array1d_iT1y_iT,depth_end):
    chlfile=os.path.join(indir,'chl.0000000000.t001.nc')
    # mg Chl-a m^-3
    array2d_idepth_iT_chlfull=\
    read_mitgcm.get_array2d_idepth_iT_chlfull(chlfile)
    # mg Chl-a m^-2
    array1d_iT_vintchlfull\
    =vstats.vint(
        array2d_idepth_iT_tracer=array2d_idepth_iT_chlfull,
        array1d_idepth_delR=drF,
        depth_end=depth_end)
    array1d_iT_vintchl=array1d_iT_vintchlfull[array1d_iT1y_iT]
    return array1d_iT_vintchl


depth_end=40
array1d_iT_modvint40mchl\
=get_array1d_iT_vintchl(indir_exp0,array1d_iT1y_iT,depth_end)


depth_end=100
array1d_iT_modvintchl\
=get_array1d_iT_vintchl(indir_exp0,array1d_iT1y_iT,depth_end)


indir_exp1='data/DataS8_output_mitgcm/exp1'
depth_end=100
array1d_iT_exp1vintchl\
=get_array1d_iT_vintchl(indir_exp1,array1d_iT1y_iT,depth_end)


indir_exp2='data/DataS8_output_mitgcm/exp2'
depth_end=100
array1d_iT_exp2vintchl\
=get_array1d_iT_vintchl(indir_exp2,array1d_iT1y_iT,depth_end)


indir_exp3times0_25='data/DataS8_output_mitgcm/exp3times0_25'
indir_exp3times0_50='data/DataS8_output_mitgcm/exp3times0_50'
indir_exp3times2_00='data/DataS8_output_mitgcm/exp3times2_00'
indir_exp3times4_00='data/DataS8_output_mitgcm/exp3times4_00'
depth_end=100
array1d_iT_exp3times0_25vintchl\
=get_array1d_iT_vintchl(indir_exp3times0_25,array1d_iT1y_iT,depth_end)
array1d_iT_exp3times0_50vintchl\
=get_array1d_iT_vintchl(indir_exp3times0_50,array1d_iT1y_iT,depth_end)
array1d_iT_exp3times2_00vintchl\
=get_array1d_iT_vintchl(indir_exp3times2_00,array1d_iT1y_iT,depth_end)
array1d_iT_exp3times4_00vintchl\
=get_array1d_iT_vintchl(indir_exp3times4_00,array1d_iT1y_iT,depth_end)


tavefile='data/DataS8_output_mitgcm/exp0/tave.0000000000.t001.nc'
array2d_idepth_iT_Ttavefull=netcdf_tools.read_netcdf(tavefile,'Ttave').squeeze().transpose()
array2d_idepth_iT_Ttavefull[-1,:]=np.nan
array2d_idepth_iT_Ttave=array2d_idepth_iT_Ttavefull[:,array1d_iT1y_iT]
minT=np.nanmin(array2d_idepth_iT_Ttave)
maxT=np.nanmax(array2d_idepth_iT_Ttave)


IPython.display.Markdown(
'''The interval of temperature ($^{\circ} C$) is
$[%f,%f].$'''
    %(minT,maxT)
)


tau_T=0.8 # unitless
A_T=-4000 # K
T_N=20+273.15 # K
mingammaT=tau_T*math.exp(A_T*(1/(minT+273.15)-1/T_N))
maxgammaT=tau_T*math.exp(A_T*(1/(maxT+273.15)-1/T_N))


IPython.display.Markdown(
'''The interval of $\gamma^T$ (unitless)
in the simulation is
$[%f,%f].$'''
    %(mingammaT,maxgammaT)
)


minpcmaxgammaT=pcmax[15-1]*mingammaT
maxpcmaxgammaT=pcmax[15-1]*maxgammaT


IPython.display.Markdown(
'''The interval of $P^C_{max}\ \gamma^T(T)$ for the diatom 6.6 $\mu$m in $d^{-1}$ is
$[%f,%f].$'''
    %(minpcmaxgammaT,maxpcmaxgammaT)
)


alphachl=info.f['alphachl'] # mmol C (mg Chl a)^-1 (umol photons m^-2)^-1
alphachl=alphachl[15-1]


modalphachl_modunits=alphachl # mmol C (mg Chl a)^-1 (umol photons m^-2)^-1
mmolC2molC=1E-3               # mol C (mmol C )^-1
molarmassC=12.011             # g C (mol C)^-1
mgChla2gChla=1E3              # mg Chl a (g Chl a)^-1
s2h=3600                      # sh^-1
modalphachl=modalphachl_modunits*mmolC2molC*molarmassC*mgChla2gChla*s2h


IPython.display.Markdown(
r'''$\alpha^{chl} = %.3G
\ g\ C\ (g\ Chl\ a)^{-1}\ h^{-1}\ (\mu mol\ photons\ m^{-2}\ s^{-1})^{-1}$.'''
    %(modalphachl)
)


array2d_idepth_iT_modPCmrep15=pcmax[15-1] \
* tau_T*np.exp(A_T*(1/(array2d_idepth_iT_Ttave+273.15)-1/T_N))


chlfile='data/DataS8_output_mitgcm/exp0/chl.0000000000.t001.nc'
array2d_idepth_iT_diatom7umchlfull\
=netcdf_tools.read_netcdf(chlfile, 'TRAC84').squeeze().transpose()
array2d_idepth_iT_diatom7umchlfull[-1,:]=np.nan
array2d_idepth_iT_diatom7umchl\
=array2d_idepth_iT_diatom7umchlfull[:,array1d_iT1y_iT]


array2d_idepth_iT_diatom7umbiofull\
=netcdf_tools.read_netcdf(carfile, 'TRAC35').squeeze().transpose()
array2d_idepth_iT_diatom7umbiofull[-1,:]=np.nan
array2d_idepth_iT_diatom7umbio\
=array2d_idepth_iT_diatom7umbiofull[:,array1d_iT1y_iT]


array2d_idepth_iT_modtheta15\
=array2d_idepth_iT_diatom7umchl/array2d_idepth_iT_diatom7umbio


s2d=86400 # s d^-1
array2d_idepth_iT_modEk15=\
array2d_idepth_iT_modPCmrep15 \
/ (modalphachl_modunits * array2d_idepth_iT_modtheta15 * s2d)


array2d_idepth_iT_modEk15light=array2d_idepth_iT_modEk15.copy()
array2d_idepth_iT_modEk15light[array2d_idepth_iT_meanPARumolEonm2s<5]=np.nan


def make_plots(axs):
    box=dict(pad=5,alpha=0)
    locs=np.array([0, 31, 59, 90, 120, 151,
              181, 212, 243, 273, 304, 334])
    labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
            'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
    
    # --- ICE CONCENTRATIONS
    ax=axs[0]
    h1=ax.plot(first_year,ice*100,'k-',lw=2.5,label='model')
    ax.set_xlim(first_year[0],first_year[-1])
    xlims=ax.get_xlim()
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_yticks([0,50,100])
    ax.yaxis.set_ticks_position('both')
    ax.set_ylabel('Ice conc.\n (%)',bbox=box)
    ax.set_ylim(0,100)
    ax.grid()
    ax.legend(bbox_to_anchor=(0.21,1.00))
    plt.text(0.03,0.7,'(a)',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=30)
    # adding the vertical line at 15% of observed ice
    ax.axvline(obs_breakup,color='k',linestyle='--')

    # --- SNOW AND ICE THICKNESS
    ax=axs[1]
    h1=ax.fill_between(first_year,array1d_iT_snow_volume,
                       linestyle='-',color='aliceblue',linewidth=1)
    h2=ax.plot(first_year,array1d_iT_snow_volume,
               linestyle='-',color='black',lw=1)
    h3=ax.plot((snow_thickness_df['doy']-1).to_numpy(),
              snow_thickness_df['sample_thickness_m_average'].to_numpy(),
              'o',color='black')
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_ylabel('Thickness\n(m)',bbox=box)
    mpl.axis.Axis.set_label_coords(ax.yaxis,-0.061,0.3)
    ax.set_ylim(-1.7,1.0)
    ax.grid()
    ax.text(60,0.2,'snow',
            horizontalalignment = 'left',verticalalignment = 'center',
            fontsize=20)
    h4=ax.fill_between(first_year,-array1d_iT_ice_volume,
                       linestyle='-',color='dodgerblue',linewidth=1)
    h5=ax.plot(first_year,-array1d_iT_ice_volume,
               linestyle='-',color='black',lw=1)
    h6=ax.plot((ice_thickness_df['doy']-1).to_numpy(),
              -ice_thickness_df['sample_thickness_m_average'].to_numpy(),
              'o',color='black')
    ax.text(60,-0.30,'ice',
            horizontalalignment = 'left',verticalalignment = 'center',
            fontsize=20)
    ax.legend([h2[0],h3[0]],['model','obs'],loc=4,
              bbox_to_anchor=(0.90,0.00))
    plt.text(0.03,0.85,'(b)',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=30)
    # adding the vertical line at 15% of observed ice
    ax.axvline(obs_breakup,linestyle='--',color='black')

    # --- sfc PAR
    ax=axs[2]
    h1=ax.plot(first_year,array1d_iT_modPAR3m,'-',color='black',lw=3)
    ax.set_yscale('log')
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_ylabel('PAR at 3 m\n(mol photons m$^{-2}$ d$^{-1})$',bbox=box)
    ax.yaxis.set_ticks_position('both')
    ax.set_ylim(1E-2,50)
    ax.grid()
    plt.text(0.03,0.90,'(c)',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=30)
    h2=ax.plot(first_year,array1d_iT_obsPAR3m,
               'o',color='black')
    ax.legend([h1[0],h2[0]],['model','obs'],loc=0, ncol=1)
    # adding the vertical line at 15% of observed ice
    ax.axvline(obs_breakup,color='k',linestyle='--')

    # --- sfc NO3 + sfc SiO2
    ax=axs[3]
    h1=ax.plot(first_year,array1d_iT_modnitrate3m,'-',color='blue',lw=3)
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_ylabel('Nutrients at 3 m\n($\mathrm{ \mu M }$)',bbox=box)
    ax.yaxis.set_ticks_position('both')
    ax.set_ylim(0,10)
    ax.grid()
    plt.text(0.03,0.90,'(d)',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=30)
    h2=ax.plot(first_year,array1d_iT_silicicacid3m ,'-',color='red', lw=3)
    h3=ax.plot(first_year,array1d_iT_obsnitrate3m,
               'o',color='blue')
    h4=ax.plot(first_year,array1d_iT_obssilica3m ,
               's',color='red')
    ax.legend([h1[0],h3[0],
               h2[0],h4[0]],
              ['model $\mathrm{ NO_3 }$','obs $\mathrm{ NO_3 }$',
               'model $\mathrm{ Si{(OH)}_4 }$','obs $\mathrm{ Si{(OH)}_4 }$'],
              loc=1, ncol=2)
    # adding the vertical line at 15% of observed ice
    ax.axvline(obs_breakup,color='k',linestyle='--')

    # --- CHLOROPHYLL VINT
    ax=axs[4]
    h1=ax.plot(first_year,array1d_iT_modvint40mchl,
               '--',lw=3,color='black',label='model (0-40m)')
    h2=ax.plot(first_year,array1d_iT_modvintchl,
               '-',color='black',lw=3,label='model (0-100m)')
    h3=ax.plot(first_year,array1d_iT_obsvint40mchlHPLC,
               '+',color='black',label='obs (0-40m)')
    h4=ax.plot(first_year,array1d_iT_obsvintchlHPLC,
               'o',color='black',label='obs (0-100m)')
    ax.set_yscale('log')
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_ylabel('Vertically integrated Chl $a$\n($\mathrm{ mg\ Chl\ m^{-2} }$)',
                  bbox=box)
    ax.yaxis.set_ticks_position('both')
    ax.set_ylim(0.2,300)
    ax.grid()
    plt.text(0.03,0.90,'(e)',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=30)
    ax.legend(bbox_to_anchor=(0.31,1.00))
    # adding the vertical line at 15% of observed ice
    ax.axvline(obs_breakup,color='k',linestyle='--')
    
    # -- ISOLUME
    ax=axs[5]
    h1=ax.plot(first_year,
               array1d_iT_modisolume,
               '-',
               lw=3,
               color='black',
               label='model')
    h2=ax.plot(first_year,
               array1d_iT_obsisolume,
               'o',
               color='black',
               label='obs')
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_ylabel('Depth of $\mathrm{ z_{0.415}\ isolume }$\n(m)',
                  bbox=box)
    ax.yaxis.set_ticks_position('both')
    ax.set_ylim(0,80)
    ax.invert_yaxis()
    ax.legend(loc=3)
    ax.grid()
    plt.text(0.03,0.90,'(f)',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=30)
    # adding the vertical line at 15% of ice
    ax.axvline(obs_breakup,color='k',linestyle='--')

    # --- ADDITIONAL X-AXIS
    ax=axs[6]
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels(labels)
    ax.axes.get_yaxis().set_visible(False)
    ax.spines['top'].set_visible(False)

    # --- POSITION
    axs[0].set_position( [0.13,0.95  ,0.85,0.03])
    axs[1].set_position( [0.13,0.88  ,0.85,0.06])
    axs[2].set_position( [0.13,0.67  ,0.85,0.20])
    axs[3].set_position( [0.13,0.46  ,0.85,0.20])
    axs[4].set_position( [0.13,0.25  ,0.85,0.20])
    axs[5].set_position( [0.13,0.04  ,0.85,0.20])
    axs[6].set_position( [0.13,0.02  ,0.85,0.01])

with plt.style.context('mplstyles/obs_and_model.mplstyle'):    
    # Plot
    fig,axs=plt.subplots(7,1,sharex=False,figsize=(16, 24))
    make_plots(axs)
    fig.align_ylabels(axs[:])

    # --- SAVE
    plt.savefig('figures_progress/qik.png',dpi=1000)


type_id=np.arange(1,52)

pos_Prochlorococcus=np.arange( 0, 1)
pos_Synechococcus  =np.arange( 1, 2)
pos_smalleuk       =np.arange( 2, 4)
pos_cocco          =np.arange( 4, 9)
pos_diazo          =np.arange( 9,13)
pos_Trichodesmium  =np.arange(13,14)
pos_diatoms        =np.arange(14,23)
pos_mixodino       =np.arange(23,33)
pos_zoo            =np.arange(33,49)


def make_plots(axs):
    box=dict(pad=5,alpha=0)
    
    # --- EQUIVALENT SPHERICAL DIAMETER (ESD)
    ax=axs[0]
    h3=ax.plot(type_id[pos_smalleuk]-2,  esd[pos_smalleuk],'gh',markersize=12)
    h4=ax.plot(type_id[pos_cocco]-2,     esd[pos_cocco],   'b+')
    h7=ax.plot(type_id[pos_diatoms]-2-5, esd[pos_diatoms], 'ro')
    h8=ax.plot(type_id[pos_mixodino]-2-5,esd[pos_mixodino],'md')
    h9=ax.plot(type_id[pos_zoo]-2-5,     esd[pos_zoo],     'k*')
    ax.set_xticklabels([])
    ax.xaxis.set_major_locator(mticker.MultipleLocator(10))
    ax.xaxis.set_minor_locator(mticker.MultipleLocator(1))
    ax.set_xlim([0,43])
    ax.set_ylabel('Equivalent spherical diameter ($\mathrm{ \mu m }$)')
    ax.set_yscale('log')
    ax.grid()
    plt.text(0.03,0.92,'(a)',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=30)

    # --- MAXIMUM GROWTH RATE (P^C_MAX)
    ax=axs[1]
    h3=ax.plot(type_id[pos_smalleuk]-2,  pcmax[pos_smalleuk],'gh',markersize=12)
    h4=ax.plot(type_id[pos_cocco]-2,     pcmax[pos_cocco],   'b+')
    h7=ax.plot(type_id[pos_diatoms]-2-5, pcmax[pos_diatoms], 'ro')
    h8=ax.plot(type_id[pos_mixodino]-2-5,pcmax[pos_mixodino],'md')
    h9=ax.plot(type_id[pos_zoo]-2-5,     pcmax[pos_zoo],     'k*')
    ax.set_xticklabels([])
    ax.xaxis.set_major_locator(mticker.MultipleLocator(10))
    ax.xaxis.set_minor_locator(mticker.MultipleLocator(1))
    ax.set_xlim([0,43])
    ax.set_ylabel('Maximum growth rate ($\mathrm{ d^{-1} }$)')
    ax.set_yscale('log')
    ax.yaxis.set_major_locator(plt.FixedLocator([0.5,1,1.5,2]))
    ax.yaxis.set_minor_locator(plt.FixedLocator([0.6,0.7,0.8,0.9,
                                                 1.1,1.2,1.3,1.4,1.6,1.7,1.8,1.9]))
    ax.yaxis.set_major_formatter(plt.FixedFormatter([0.5,1,1.5,2]))
    ax.yaxis.set_minor_formatter(plt.FixedFormatter(['']))
    ax.set_ylim([0.5,2])
    ax.grid()
    plt.text(0.03,0.92,'(b)',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=30)

    # --- HALF SATURATION FOR GROWTH ON NITRATE (k_NO_3)
    ax=axs[2]
    h3=ax.plot(type_id[pos_smalleuk]-2,  ksatno3[pos_smalleuk],'gh',markersize=12)
    h4=ax.plot(type_id[pos_cocco]-2,     ksatno3[pos_cocco],   'b+')
    h7=ax.plot(type_id[pos_diatoms]-2-5, ksatno3[pos_diatoms], 'ro')
    h8=ax.plot(type_id[pos_mixodino]-2-5,ksatno3[pos_mixodino],'md')
    h9=ax.plot(type_id[pos_zoo]-2-5,     ksatno3[pos_zoo],     'k*')
    ax.set_xticklabels([])
    ax.xaxis.set_major_locator(mticker.MultipleLocator(10))
    ax.xaxis.set_minor_locator(mticker.MultipleLocator(1))
    ax.set_xlim([0,43])
    ax.set_ylabel('$\mathrm{ k_{NO_3} }$ ($\mathrm{ mmol\ N\ m^{-3} }$)')
    ax.set_yscale('log')
    ax.set_ylim([0.01,10])
    ax.grid()
    plt.text(0.03,0.92,'(c)',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=30)

    # --- MAXIMUM GRAZING RATE (G_MAX)
    ax=axs[3]
    h3=ax.plot(type_id[pos_smalleuk]-2,  gmax[pos_smalleuk],'gh',markersize=12,
               label='picophytoplankton')
    h4=ax.plot(type_id[pos_cocco]-2,     gmax[pos_cocco],   'b+',
               label='other nanophytoplankton')
    h7=ax.plot(type_id[pos_diatoms]-2-5, gmax[pos_diatoms], 'ro',
               label='diatoms')
    h8=ax.plot(type_id[pos_mixodino]-2-5,gmax[pos_mixodino],'md',
               label='mixotrophic dinoflagellates')
    h9=ax.plot(type_id[pos_zoo]-2-5,     gmax[pos_zoo],     'k*',
               label='microzooplankton')
    ax.set_xlabel('Planktonic type ID',fontsize=23)
    ax.xaxis.set_major_locator(mticker.MultipleLocator(10))
    ax.xaxis.set_minor_locator(mticker.MultipleLocator(1))
    ax.set_xlim([0,43])
    ax.set_ylabel('Maximum grazing rate ($\mathrm{ d^{-1} }$)')
    ax.set_yscale('log')
    ax.yaxis.set_major_locator(plt.FixedLocator([0.5,1,5,10]))
    ax.yaxis.set_minor_locator(plt.FixedLocator([0.6,0.7,0.8,0.9,2,3,4,6,7,8,9,20]))
    ax.yaxis.set_major_formatter(plt.FixedFormatter([0.5,1,5,10]))
    ax.set_ylim([0.5,20])
    ax.grid()
    ax.legend(loc=3)
    plt.text(0.03,0.92,'(d)',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=30)

    # --- POSITION
    axs[0].set_position( [0.10,0.76  ,0.85,0.21])
    axs[1].set_position( [0.10,0.53  ,0.85,0.21])
    axs[2].set_position( [0.10,0.30  ,0.85,0.21])
    axs[3].set_position( [0.10,0.07  ,0.85,0.21])

with plt.style.context('mplstyles/gud_traits.mplstyle'):    
    # Plot
    fig,axs=plt.subplots(4,1,figsize=(16, 16))
    make_plots(axs)
    fig.align_ylabels(axs[:])

    # --- SAVE
    plt.savefig('figures_progress/traits.png',dpi=1000)


def make_plots(axs):
    locs=np.array([0, 31, 59, 90, 120, 151,
              181, 212, 243, 273, 304, 334])
    labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
            'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
    
    # --- NITRATE
    ax=axs[0]
    h=ax.pcolormesh(first_year366,
                    RF_above151,
                    array2d_idepth_iT_modno3[0:(RF_above151.size)-1,:],
                    cmap='viridis',
                    vmin=0,
                    vmax=6)
    ax.scatter(nutrients_df['doy']-1,
               nutrients_df['depth'],
               c=nutrients_df['no3'],
               cmap='viridis',
               vmin=0,
               vmax=6,
               s=100,
               edgecolor='k')
    ax.set_xticklabels([])
    ax.set_ylabel('Depth (m)')
    ax.set_xlim(first_year366[0],first_year[-1])
    ax.set_ylim(150,0)
    ax.grid()
    ax.set_xticks(locs)
    ax.set_xticklabels(labels)
    cbaxes = fig.add_axes([0.88, 0.83, 0.01, 0.14])
    cbar=plt.colorbar(h,
                      cax = cbaxes)
    cbar.set_label('(mmol N $\mathrm{ m^{-3} }$)')
    plt.text(0.50,1.10,'(a) nitrate',
             horizontalalignment = 'center', verticalalignment='center',
             transform=ax.transAxes)
    
    # --- SILICIC ACID
    ax=axs[1]
    h=ax.pcolormesh(first_year366,
                    RF_above151,
                    array2d_idepth_iT_modsioh4[0:(RF_above151.size)-1,:],
                    cmap='viridis',
                    vmin=0,
                    vmax=10)
    ax.scatter(nutrients_df['doy']-1,
               nutrients_df['depth'],
               c=nutrients_df['sioh4'],
               cmap='viridis',
               vmin=0,
               vmax=10,
               s=100,
               edgecolor='k')
    ax.set_xticklabels([])
    ax.set_ylabel('Depth (m)')
    ax.set_xlim(first_year366[0],first_year[-1])
    ax.set_ylim(150,0)
    ax.grid()
    ax.set_xticks(locs)
    ax.set_xticklabels(labels)
    cbaxes = fig.add_axes([0.88, 0.63, 0.01, 0.14])
    cbar=plt.colorbar(h,
                      cax = cbaxes)
    cbar.set_label('(mmol Si$\mathrm{ {(OH)}_4 }$ $\mathrm{ m^{-3} }$)')
    plt.text(0.50,1.10,'(b) silicic acid',
             horizontalalignment = 'center', verticalalignment='center',
             transform=ax.transAxes)
    
    # --- PHOSPHATE
    ax=axs[2]
    h=ax.pcolormesh(first_year366,
                    RF_above151,
                    array2d_idepth_iT_modpo4[0:(RF_above151.size)-1,:],
                    cmap='viridis',
                    vmin=0,
                    vmax=1)
    ax.scatter(nutrients_df['doy']-1,
               nutrients_df['depth'],
               c=nutrients_df['po4'],
               cmap='viridis',
               vmin=0,
               vmax=1,
               s=100,
               edgecolor='k')
    ax.set_xticklabels([])
    ax.set_ylabel('Depth (m)')
    ax.set_xlim(first_year366[0],first_year[-1])
    ax.set_ylim(150,0)
    ax.grid()
    ax.set_xticks(locs)
    ax.set_xticklabels(labels)
    cbaxes = fig.add_axes([0.88, 0.43, 0.01, 0.14])
    cbar=plt.colorbar(h,
                      cax = cbaxes)
    cbar.set_label('(mmol $\mathrm{ PO_4 }$ $\mathrm{ m^{-3} }$)')
    plt.text(0.50,1.10,'(c) phosphate',
             horizontalalignment = 'center', verticalalignment='center',
             transform=ax.transAxes)
    
    # --- IRON
    ax=axs[3]
    h=ax.pcolormesh(first_year366,
                    RF_above151,
                    array2d_idepth_iT_modfet[0:(RF_above151.size)-1,:],
                    cmap='viridis',
                    vmin=0,
                    vmax=1.5E-3)
    ax.scatter(fet_df['doy']-1,
               fet_df['depth'],
               c=fet_df['Fe_D_CONC_BOTTLE_mmol_m3'],
               cmap='viridis',
               vmin=0,
               vmax=1.5E-3,
               s=100,
               edgecolor='k')
    ax.set_xticklabels([])
    ax.set_ylabel('Depth (m)')
    ax.set_xlim(first_year366[0],first_year[-1])
    ax.set_ylim(150,0)
    ax.grid()
    ax.set_xticks(locs)
    ax.set_xticklabels(labels)
    cbaxes = fig.add_axes([0.88, 0.23, 0.01, 0.14])
    cbar=plt.colorbar(h,
                      cax = cbaxes)
    cbar.set_label('(mmol Fe $\mathrm{ m^{-3} }$)')
    plt.text(0.50,1.10,'(d) total dissolved inorganic iron',
             horizontalalignment = 'center', verticalalignment='center',
             transform=ax.transAxes)
    
    # --- CHLOROPHYLL A
    ax=axs[4]
    h=ax.pcolormesh(first_year366,
                    RF_above151,
                    array2d_idepth_iT_modchl[0:(RF_above151.size)-1,:],
                    cmap='viridis',
                    norm=mpl.colors.LogNorm(vmin=1E-2, vmax=10))
    ax.scatter(chlHPLC_df['doy']-1,
               chlHPLC_df['depth'],
               c=chlHPLC_df['chlHPLC'],
               cmap='viridis',
               norm=mpl.colors.LogNorm(vmin=1E-2, vmax=10),
               s=100,
               edgecolor='k')
    ax.set_xticklabels([])
    ax.set_ylabel('Depth (m)')
    ax.set_xlim(first_year366[0],first_year[-1])
    ax.set_ylim(150,0)
    ax.grid()
    ax.set_xticks(locs)
    ax.set_xticklabels(labels)
    cbaxes = fig.add_axes([0.88, 0.03, 0.01, 0.14])
    cbar=plt.colorbar(h,
                      cax = cbaxes)
    cbar.set_label('(mg Chl $a$ $\mathrm{ m^{-3} }$)')
    plt.text(0.50,1.10,'(e) chlorophyll $a$',
             horizontalalignment = 'center', verticalalignment='center',
             transform=ax.transAxes)

    # --- POSITION
    axs[0].set_position( [0.07,0.83  ,0.79,0.14])
    axs[1].set_position( [0.07,0.63  ,0.79,0.14])
    axs[2].set_position( [0.07,0.43  ,0.79,0.14])
    axs[3].set_position( [0.07,0.23  ,0.79,0.14])
    axs[4].set_position( [0.07,0.03  ,0.79,0.14])

with plt.style.context('mplstyles/heatmaps.mplstyle'):    
    # Plot
    fig,axs=plt.subplots(5,1,sharex=False,figsize=(16, 24))
    make_plots(axs)

    # --- SAVE
    plt.savefig('figures_progress/conc.png',dpi=1000)


obsPCmax_df=obsPCmax_df.copy() # to avoid a SettingWithCopyWarning
obsPCmax_df.sort_values(by=['species'],inplace=True)
obsalphachl_df=obsalphachl_df.copy() # to avoid a SettingWithCopyWarning
obsalphachl_df.sort_values(by=['species'],inplace=True)
obsEk_df=obsEk_df.copy() # to avoid a SettingWithCopyWarning
obsEk_df.sort_values(by=['species'],inplace=True)
array1d_i_modEk15light=array2d_idepth_iT_modEk15light.flatten()
pos=~np.isnan(array1d_i_modEk15light)
array1d_i_modEk15light=array1d_i_modEk15light[pos]


def make_plots(axs):
    
    # --- MAXIMUM PHOTOSYNTHESIS RATE AT SATURATION IRRADIANCE
    # --- IN NUTRIENT REPLETE CONDITIONS
    # --- (P^C_max * gamma^T)
    ax=axs[0,0]
    xlim=(0,1.5)
    h1=ax.scatter(obsPCmax_df.value,obsPCmax_df.species,
                  color='b',label='obs')
    h2=ax.axhline('',
            xmin=minpcmaxgammaT/xlim[1],
            xmax=maxpcmaxgammaT/xlim[1],
            color='red',lw=10)
    h3=ax.axhline('foo',
                 lw=0)
    # to have the ytick labels in alphabetical order from top to bottom
    ax.invert_yaxis()
    # in italics
    # FULL GENUS
    ax.set_yticklabels(['Chaetoceros furcellatus',
                        'Chaetoceros neogracile',
                        'Fragilariopsis cylindrus',
                        'Thalassiosira gravida 1',
                        'Thalassiosira gravida 2',
                        'Thalassiosira nordenskioeldii',
                        'Numerical diatom 6.6 μm'
                       ],
                       fontstyle='italic',
                       fontweight='bold')
    ax.set_xlim(xlim)
    ax.set_xlabel('$P^C_{max} \gamma^T$\n($\mathrm{ d^{-1} }$)')
    ax.grid()
    plt.text(0.10,0.88,'(a)',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=30)
    ax.legend([h1,h2],['obs','model'],ncol=2,loc='lower center')
    
    # --- NOTHING
    axs[1,0].axis('off')
    
    # --- LINEAR INITIAL SLOPE OF THE CHL A-SPECIFIC PHOTOSYNTHESIS
    # --- VERSUS IRRADIANCE CURVE
    # --- (ALPHA^CHL)
    ax=axs[0,1]
    xlim=(0,0.05)
    h1=ax.scatter(obsalphachl_df.value,
                  obsalphachl_df.species,
                  c='b',label='obs')
    h2=ax.scatter(modalphachl,
                  '',
                  color='red')
    h3=ax.axhline('foo',
                  lw=0)
    # to have the ytick labels in alphabetical order from top to bottom
    ax.invert_yaxis()
    ax.set_yticklabels([])
    ax.set_xlim(xlim)
    ax.set_xlabel('''$\\alpha^{chl}$\n\
    $\mathrm{ (g\ C\ (g\ Chl)^{-1}\ h^{-1} }$\n\
    $\mathrm{ (\mu mol\ photons\ m^{-2}\ s^{-1})^{-1}) }$''')
    ax.grid()
    plt.text(0.10,0.88,'(b)',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=30)
    
    # --- NOTHING
    axs[1,1].axis('off')
    
    # --- LIGHT SATURATION PARAMETER
    # --- (E_k)
    ax=axs[0,2]
    xlim=(0,250)
    h1=ax.scatter(obsEk_df.value,
                  obsEk_df.species,
                  c='b',label='obs')
    h3=ax.axhline('foo',
                  lw=0)
    h4=ax.axhline('foo2',
                  lw=0)
    # to have the ytick labels in alphabetical order from top to bottom
    ax.invert_yaxis()
    ax.set_yticklabels([])
    ax.set_xlabel('$E_k$\n$\mathrm{ (\mu mol\ photons\ m^{-2}\ s^{-1}) }$')
    ax.set_xlim(xlim)
    ax.grid()
    plt.text(0.10,0.88,'(c)',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=30)

    ax=axs[1,2]
    boxplot=ax.boxplot(array1d_i_modEk15light,
                       vert=False,
                       widths=0.75,
                       boxprops=dict(color='red'),
                       flierprops=dict(marker='o',
                                       markerfacecolor='red',
                                       markeredgecolor='none',
                                       markersize=1),
                       medianprops=dict(color='red'),
                       capprops=dict(color='red'),
                       whiskerprops=dict(color='red'))
    ax.set_xlim(xlim)
    ax.set_yticklabels([])
    ax.set_xticklabels([])
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_visible(False)
#     ax.spines['left'].set_visible(False)
    ax.tick_params(length=0)
    ax.grid()
    
    # --- POSITION
    # FULL GENUS
    axs[0,0].set_position( [0.24,0.21  ,0.22,0.77])
    axs[0,1].set_position( [0.50,0.21  ,0.22,0.77])
    axs[0,2].set_position( [0.76,0.21  ,0.22,0.77])
    axs[1,2].set_position( [0.76,0.31  ,0.22,0.07])
    
    return boxplot

with plt.style.context('mplstyles/lacour17.mplstyle'):    
    # Plot
    fig,axs=plt.subplots(2,3,figsize=(16, 8))
    boxplot=make_plots(axs)

    # --- SAVE
    plt.savefig('figures_progress/lacour17.png',dpi=1000)


upper_whisker=[item.get_xdata()[1] for item in boxplot['whiskers']][1]
fliers=array1d_i_modEk15light[array1d_i_modEk15light>upper_whisker]
nb_fliers=fliers.size
nb_total=array1d_i_modEk15light.size


nb_fliers

311


nb_total

3409


data={'exp0x':[0,    163,  198,228,300,328,364],
      'exp0y':[1,    1,    100,100, 10,  1,  1],
      'exp1x':[0,    201,  225,255,300,328,364],
      'exp1y':[1,    1,    100,100, 10,  1,  1],
      'exp2x':[0,    163,  228,258,300,328,364],
      'exp2y':[0.001,0.001,100,100, 10,  1,  1],
      'exp3x':[0,    163,  181,210,300,325,364],
      'exp3y':[1,    1,     10, 10,  3,  1,  1]}
df=pd.DataFrame(data)
# df


def make_plots(axs):
    locs=np.array([0, 31, 59, 90, 120, 151,
              181, 212, 243, 273, 304, 334])
    labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
            'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
    xlims=(first_year[0],first_year[-1])
    xmin=xlims[0]
    xmax=xlims[1]
    
    # --- EXP-0: ICE CONCENTRATION
    ax=axs[0,0]
    h1=ax.plot(first_year,ice*100,'k-',lw=2.5,label='model')
    ax.set_xlim(xlims)
    ymin=0
    ymax=100
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_yticks([0,50,100])
    ax.yaxis.set_ticks_position('both')
    ax.set_ylabel('Ice conc.\n (%)')
    ax.set_ylim(0,100)
    ax.grid()
    ax.legend(loc='lower left',bbox_to_anchor=(0.08,-0.15))
    ax.axvline(obs_breakup,color='k',linestyle='--')
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    plt.text(0.03,0.60,'(a)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    plt.text(0.5,1.5,'EXP-0: default',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    # --- EXP-1: ICE CONCENTRATION
    ax=axs[0,1]
    h1=ax.plot(first_year,ice*100,'k-',lw=2.5,label='model')
    ax.set_xlim(xlims)
    ymin=0
    ymax=100
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_yticks([0,50,100])
    ax.yaxis.set_ticks_position('both')
    ax.set_ylabel('Ice conc.\n (%)')
    ax.set_ylim(ymin,ymax)
    ax.grid()
    ax.legend(loc='lower left',bbox_to_anchor=(0.08,-0.15))
    ax.axvline(obs_breakup,color='k',linestyle='--')
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    plt.text(0.03,0.60,'(b)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    plt.text(0.5,1.5,'EXP-1: no light under sea ice',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)

    # --- EXP-0: CHLOROPHYLL VINT
    ax=axs[1,0]
    h1=ax.plot(first_year,array1d_iT_modvintchl,
               '-',color='black',lw=3,label='model (0-100m)')
    h2=ax.plot(first_year,array1d_iT_obsvintchlHPLC,
               'o',color='black',label='obs (0-100m)')
    ax.set_yscale('log')
    ymin=0.2
    ymax=300
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_ylabel('Vertically integrated Chl $a$\n($\mathrm{ mg\ Chl\ m^{-2} }$)')
    ax.yaxis.set_ticks_position('both')
    ax.set_ylim(ymin,ymax)
    ax.grid()
    ax.legend(loc='lower left',bbox_to_anchor=(0.08,0.75))
    ax.axvline(obs_breakup,color='k',linestyle='--')
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    plt.text(0.03,0.90,'(c)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)

    # --- EXP-1: CHLOROPHYLL VINT
    ax=axs[1,1]
    h1=ax.plot(first_year,array1d_iT_exp1vintchl,
               '-.',color='red',lw=3,label='model (0-100m)')
    h2=ax.plot(first_year,array1d_iT_obsvintchlHPLC,
               'o',color='black',label='obs (0-100m)')
    ax.set_yscale('log')
    ymin=0.2
    ymax=300
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_ylabel('Vertically integrated Chl $a$\n($\mathrm{ mg\ Chl\ m^{-2} }$)')
    ax.yaxis.set_ticks_position('both')
    ax.set_ylim(ymin,ymax)
    ax.grid()
    ax.legend(loc='lower left',bbox_to_anchor=(0.08,0.75))
    ax.axvline(obs_breakup,color='k',linestyle='--')
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    plt.text(0.03,0.90,'(d)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    # --- EXP-0: ADDITIONAL X-AXIS
    ax=axs[2,0]
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels(labels)
    ax.axes.get_yaxis().set_visible(False)
    ax.spines['top'].set_visible(False)
    
    # --- EXP-1: ADDITIONAL X-AXIS
    ax=axs[2,1]
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels(labels)
    ax.axes.get_yaxis().set_visible(False)
    ax.spines['top'].set_visible(False)
    
    # --- EXP-0: SCHEMATIC
    ax=axs[3,0]
    h0=ax.plot(np.array(df.exp0x),np.array(df.exp0y),
               '-',color='black',lw=3,label='EXP-0')
    ax.set_xlabel('Time')
    ax.xaxis.set_label_coords(0.05,-0.05)
    ax.set_xlim(xlims)
    xticks=([df.exp0x.iloc[1],df.exp0x.iloc[2]])
    xticklabels=[r'$\leftarrow t_I \rightarrow$',
                r'$\leftarrow t_P \rightarrow$']
    ax.set_xticks(xticks)
    ax.set_xticklabels(xticklabels)
    ax.set_yscale('log')
    ax.set_ylabel('Log10(Chlorophyll $a$)')
    ax.yaxis.set_label_coords(-0.075,0.50)
    ymin=1E-4
    ymax=300
    ax.set_ylim(ymin,ymax)
    ax.tick_params(axis='y',which='both',left=False,labelleft=False)
    ax.plot(df.exp0x.iloc[1],df.exp0y.iloc[1],
           'o',color='black')
    ax.text(df.exp0x.iloc[1]-6,df.exp0y.iloc[1]+1,'I')
    ax.plot(df.exp0x.iloc[2],df.exp0y.iloc[2],
           'o',color='black')
    ax.text(df.exp0x.iloc[2]-6,df.exp0y.iloc[2]+20,'P')
    ax.axvline(df.exp0x.iloc[1],
               ymax=( math.log10(df.exp0y.iloc[1])-math.log10(ymin) ) \
               /    ( math.log10(ymax)            -math.log10(ymin) ),
               color='black',linestyle='--')
    ax.axvline(df.exp0x.iloc[2],
               ymax=( math.log10(df.exp0y.iloc[2])-math.log10(ymin) ) \
               /    ( math.log10(ymax)            -math.log10(ymin) ),
               color='black',linestyle='--')
    
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    
    ax.legend(loc='upper right')
    plt.text(0.03,0.90,'(e)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    # --- EXP-1: SCHEMATIC
    ax=axs[3,1]
    h0=ax.plot(np.array(df.exp1x),np.array(df.exp1y),
               '-.',color='red',lw=3,label='EXP-1')
    ax.set_xlabel('Time')
    ax.xaxis.set_label_coords(0.05,-0.05)
    ax.set_xlim(xlims)
    xticks=([df.exp1x.iloc[1],df.exp1x.iloc[2]])
    xticklabels=[r'$\leftarrow t_I \rightarrow$',
                r'$\leftarrow t_P \rightarrow$']
    ax.set_xticks(xticks)
    ax.set_xticklabels(xticklabels)
    ax.set_yscale('log')
    ax.set_ylabel('Log10(Chlorophyll $a$)')
    ax.yaxis.set_label_coords(-0.075,0.50)
    ymin=1E-4
    ymax=300
    ax.set_ylim(ymin,ymax)
    ax.tick_params(axis='y',which='both',left=False,labelleft=False)
    ax.plot(df.exp1x.iloc[1],df.exp1y.iloc[1],
           'o',color='red')
    ax.text(df.exp1x.iloc[1]-6,df.exp1y.iloc[1]+1,'I\'',color='red')
    ax.plot(df.exp1x.iloc[2],df.exp1y.iloc[2],
           'o',color='red')
    ax.text(df.exp1x.iloc[2]-6,df.exp1y.iloc[2]+20,'P\'',color='red')
    ax.axvline(df.exp1x.iloc[1],
               ymax=( math.log10(df.exp1y.iloc[1])-math.log10(ymin) ) \
               /    ( math.log10(ymax)            -math.log10(ymin) ),
               color='black',linestyle='--')
    ax.axvline(df.exp1x.iloc[2],
               ymax=( math.log10(df.exp1y.iloc[2])-math.log10(ymin) ) \
               /    ( math.log10(ymax)            -math.log10(ymin) ),
               color='black',linestyle='--')
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    
    ax.legend(loc='upper right')
    plt.text(0.03,0.90,'(f)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    # --- POSITION
    axs[0,0].set_position( [0.05,0.90  ,0.43,0.05])
    axs[0,1].set_position( [0.54,0.90  ,0.43,0.05])
    axs[1,0].set_position( [0.05,0.49  ,0.43,0.38])
    axs[1,1].set_position( [0.54,0.49  ,0.43,0.38])
    axs[2,0].set_position( [0.05,0.47  ,0.43,0.01])
    axs[2,1].set_position( [0.54,0.47  ,0.43,0.01])
    axs[3,0].set_position( [0.05,0.04  ,0.43,0.38])
    axs[3,1].set_position( [0.54,0.04  ,0.43,0.38])

with plt.style.context('mplstyles/exp.mplstyle'):    
    # Plot
    fig,axs=plt.subplots(4,2,figsize=(31.6,12.3))
    make_plots(axs)

    # --- SAVE
    plt.savefig('figures_progress/exp1.png',dpi=1000)


def make_plots(axs):
    locs=np.array([0, 31, 59, 90, 120, 151,
              181, 212, 243, 273, 304, 334])
    labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
            'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
    xlims=(first_year[0],first_year[-1])
    xmin=xlims[0]
    xmax=xlims[1]
    
    # --- EXP-0: ICE CONCENTRATION
    ax=axs[0,0]
    h1=ax.plot(first_year,ice*100,'k-',lw=2.5,label='model')
    ax.set_xlim(xlims)
    ymin=0
    ymax=100
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_yticks([0,50,100])
    ax.yaxis.set_ticks_position('both')
    ax.set_ylabel('Ice conc.\n (%)')
    ax.set_ylim(0,100)
    ax.grid()
    ax.legend(loc='lower left',bbox_to_anchor=(0.08,-0.15))
    ax.axvline(obs_breakup,color='k',linestyle='--')
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    plt.text(0.03,0.60,'(a)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    plt.text(0.5,1.5,'EXP-0: default',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    # --- EXP-2: ICE CONCENTRATION
    ax=axs[0,1]
    h1=ax.plot(first_year,ice*100,'k-',lw=2.5,label='model')
    ax.set_xlim(xlims)
    ymin=0
    ymax=100
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_yticks([0,50,100])
    ax.yaxis.set_ticks_position('both')
    ax.set_ylabel('Ice conc.\n (%)')
    ax.set_ylim(ymin,ymax)
    ax.grid()
    ax.legend(loc='lower left',bbox_to_anchor=(0.08,-0.15))
    ax.axvline(obs_breakup,color='k',linestyle='--')
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    plt.text(0.03,0.60,'(b)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    plt.text(0.5,1.5,'EXP-2: no minimum biomass',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)

    # --- EXP-0: CHLOROPHYLL VINT
    ax=axs[1,0]
    h1=ax.plot(first_year,array1d_iT_modvintchl,
               '-',color='black',lw=3,label='model (0-100m)')
    h2=ax.plot(first_year,array1d_iT_obsvintchlHPLC,
               'o',color='black',label='obs (0-100m)')
    ax.set_yscale('log')
    ymin=1E-4
    ymax=300
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_ylabel('Vertically integrated Chl $a$\n($\mathrm{ mg\ Chl\ m^{-2} }$)')
    ax.yaxis.set_ticks_position('both')
    ax.set_ylim(ymin,ymax)
    ax.grid()
    ax.legend(loc='lower left',bbox_to_anchor=(0.08,0.75))
    ax.axvline(obs_breakup,color='k',linestyle='--')
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    plt.text(0.03,0.90,'(c)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)

    # --- EXP-2: CHLOROPHYLL VINT
    ax=axs[1,1]
    h1=ax.plot(first_year,array1d_iT_exp2vintchl,
               '--',color='blue',lw=3,label='model (0-100m)')
    h2=ax.plot(first_year,array1d_iT_obsvintchlHPLC,
               'o',color='black',label='obs (0-100m)')
    ax.set_yscale('log')
    ymin=1E-4
    ymax=300
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_ylabel('Vertically integrated Chl $a$\n($\mathrm{ mg\ Chl\ m^{-2} }$)')
    ax.yaxis.set_ticks_position('both')
    ax.set_ylim(ymin,ymax)
    ax.grid()
    ax.legend(loc='lower left',bbox_to_anchor=(0.08,0.75))
    ax.axvline(obs_breakup,color='k',linestyle='--')
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    plt.text(0.03,0.90,'(d)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    # --- EXP-0: ADDITIONAL X-AXIS
    ax=axs[2,0]
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels(labels)
    ax.axes.get_yaxis().set_visible(False)
    ax.spines['top'].set_visible(False)
    
    # --- EXP-2: ADDITIONAL X-AXIS
    ax=axs[2,1]
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels(labels)
    ax.axes.get_yaxis().set_visible(False)
    ax.spines['top'].set_visible(False)
    
    # --- EXP-0: SCHEMATIC
    ax=axs[3,0]
    h0=ax.plot(np.array(df.exp0x),np.array(df.exp0y),
               '-',color='black',lw=3,label='EXP-0')
    ax.set_xlabel('Time')
    ax.xaxis.set_label_coords(0.05,-0.05)
    ax.set_xlim(xlims)
    xticks=([df.exp0x.iloc[2]])
    xticklabels=[r'$\leftarrow t_P \rightarrow$']
    ax.set_xticks(xticks)
    ax.set_xticklabels(xticklabels)
    ax.set_yscale('log')
    ax.set_ylabel('Log10(Chlorophyll $a$)')
    ax.yaxis.set_label_coords(-0.075,0.50)
    ymin=1E-4
    ymax=300
    ax.set_ylim(ymin,ymax)
    ax.tick_params(axis='y',which='minor',left=False,labelleft=False)
    yticks=([df.exp0y.iloc[0]])
    yticklabels=[r'$\uparrow$'
                 '\n'
                 r'$y_I$'
                '\n'
                r'$\downarrow$']
    ax.set_yticks(yticks)
    ax.set_yticklabels(yticklabels)
    ax.plot(df.exp0x.iloc[1],df.exp0y.iloc[1],
           'o',color='black')
    ax.text(df.exp0x.iloc[1]-6,df.exp0y.iloc[1]+1,'I')
    ax.plot(df.exp0x.iloc[2],df.exp0y.iloc[2],
           'o',color='black')
    ax.text(df.exp0x.iloc[2]-6,df.exp0y.iloc[2]+20,'P')
    ax.axvline(df.exp0x.iloc[2],
               ymax=( math.log10(df.exp0y.iloc[2])-math.log10(ymin) ) \
               /    ( math.log10(ymax)            -math.log10(ymin) ),
               color='black',linestyle='--')
    
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    
    ax.legend(loc='upper right')
    plt.text(0.03,0.90,'(e)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    # --- EXP-2: SCHEMATIC
    ax=axs[3,1]
    h0=ax.plot(np.array(df.exp2x),np.array(df.exp2y),
               '--',color='blue',lw=3,label='EXP-1')
    ax.set_xlabel('Time')
    ax.xaxis.set_label_coords(0.05,-0.05)
    ax.set_xlim(xlims)
    xticks=([df.exp2x.iloc[2]])
    xticklabels=[r'$\leftarrow t_P \rightarrow$']
    ax.set_xticks(xticks)
    ax.set_xticklabels(xticklabels)
    ax.set_yscale('log')
    ax.set_ylabel('Log10(Chlorophyll $a$)')
    ax.yaxis.set_label_coords(-0.075,0.50)
    ymin=1E-4
    ymax=300
    ax.set_ylim(ymin,ymax)
    ax.tick_params(axis='y',which='minor',left=False,labelleft=False)
    yticks=([df.exp2y.iloc[0]])
    yticklabels=[r'$\uparrow$'
                 '\n'
                 r'$y_I$'
                '\n'
                r'$\downarrow$']
    ax.set_yticks(yticks)
    ax.set_yticklabels(yticklabels)
    ax.plot(df.exp2x.iloc[1],df.exp2y.iloc[1],
           'o',color='blue')
    ax.text(df.exp2x.iloc[1]-6,df.exp2y.iloc[1]+1E-3,'I\'\'',color='blue')
    ax.plot(df.exp2x.iloc[2],df.exp2y.iloc[2],
           'o',color='blue')
    ax.text(df.exp2x.iloc[2]-6,df.exp2y.iloc[2]+20,'P\'\'',color='blue')
    ax.axvline(df.exp2x.iloc[2],
               ymax=( math.log10(df.exp2y.iloc[2])-math.log10(ymin) ) \
               /    ( math.log10(ymax)            -math.log10(ymin) ),
               color='black',linestyle='--')
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    
    ax.legend(loc='upper right')
    plt.text(0.03,0.90,'(f)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    # --- POSITION
    axs[0,0].set_position( [0.05,0.90  ,0.43,0.05])
    axs[0,1].set_position( [0.54,0.90  ,0.43,0.05])
    axs[1,0].set_position( [0.05,0.49  ,0.43,0.38])
    axs[1,1].set_position( [0.54,0.49  ,0.43,0.38])
    axs[2,0].set_position( [0.05,0.47  ,0.43,0.01])
    axs[2,1].set_position( [0.54,0.47  ,0.43,0.01])
    axs[3,0].set_position( [0.05,0.04  ,0.43,0.38])
    axs[3,1].set_position( [0.54,0.04  ,0.43,0.38])

with plt.style.context('mplstyles/exp.mplstyle'):    
    # Plot
    fig,axs=plt.subplots(4,2,figsize=(31.6,12.3))
    make_plots(axs)

    # --- SAVE
    plt.savefig('figures_progress/exp2.png',dpi=1000)


def make_plots(axs):
    locs=np.array([0, 31, 59, 90, 120, 151,
              181, 212, 243, 273, 304, 334])
    labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
            'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
    xlims=(first_year[0],first_year[-1])
    xmin=xlims[0]
    xmax=xlims[1]
    
    # --- EXP-0: ICE CONCENTRATION
    ax=axs[0]
    h1=ax.plot(first_year,ice*100,'k-',lw=2.5,label='model')
    ax.set_xlim(xlims)
    ymin=0
    ymax=100
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_yticks([0,50,100])
    ax.yaxis.set_ticks_position('both')
    ax.set_ylabel('Ice conc.\n (%)')
    ax.set_ylim(0,100)
    ax.grid()
    ax.legend(loc='lower left',bbox_to_anchor=(0.08,-0.15))
    ax.axvline(obs_breakup,color='k',linestyle='--')
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    plt.text(0.03,0.60,'(a)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    plt.text(0.5,1.5,'EXP-3: differing wintertime nitrate',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)

    # --- EXP-0: CHLOROPHYLL VINT
    ax=axs[1]
    h0_25=ax.plot(first_year,array1d_iT_exp3times0_25vintchl,
                  ':',color='violet',lw=3,label='* 1/4')
    h0_50=ax.plot(first_year,array1d_iT_exp3times0_50vintchl,
                  '-',color='blue',lw=3,label='* 1/2')
    h1_00=ax.plot(first_year,array1d_iT_modvintchl,
                  '-',color='black',lw=3,label='default')
    h2_00=ax.plot(first_year,array1d_iT_exp3times2_00vintchl,
                  '-',color='orange',lw=3,label='* 2')
    h4_00=ax.plot(first_year,array1d_iT_exp3times4_00vintchl,
                  ':',color='red',lw=3,label='${[NO_3]}^{winter}$ * 4')
    ax.set_yscale('log')
    ymin=1E0
    ymax=1E2
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels([])
    ax.set_ylabel('Vertically integrated Chl $a$\n($\mathrm{ mg\ Chl\ m^{-2} }$)')
    ax.yaxis.set_ticks_position('both')
    ax.set_ylim(ymin,ymax)
    ax.grid()
    #get handles and labels
    handles, legend_labels = ax.get_legend_handles_labels()
    #specify order of items in legend
    order = [2,4,3,1,0]
    #add legend to plot
    ax.legend([handles[idx] for idx in order],
              [legend_labels [idx] for idx in order],
               loc='upper left',
               bbox_to_anchor=(0.08,1.00))
    ax.axvline(obs_breakup,color='k',linestyle='--')
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    plt.text(0.03,0.90,'(b)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    # --- EXP-0: ADDITIONAL X-AXIS
    ax=axs[2]
    ax.set_xlim(xlims)
    ax.set_xticks(locs)
    ax.set_xticklabels(labels)
    ax.axes.get_yaxis().set_visible(False)
    ax.spines['top'].set_visible(False)
    
    # --- EXP-0: SCHEMATIC
    ax=axs[3]
    h0=ax.plot(np.array(df.exp0x),np.array(df.exp0y),
               '-',color='black',lw=3,label='EXP-0')
    ax.set_xlabel('Time')
    ax.xaxis.set_label_coords(0.05,-0.05)
    ax.set_xlim(xlims)
    xticks=([df.exp0x.iloc[2]])
    xticklabels=[r'$\leftarrow t_P \rightarrow$']
    ax.set_xticks(xticks)
    ax.set_xticklabels(xticklabels)
    ax.set_yscale('log')
    ax.set_ylabel('Log10(Chlorophyll $a$)')
    ax.yaxis.set_label_coords(-0.075,0.50)
    ymin=1E-4
    ymax=300
    ax.set_ylim(ymin,ymax)
    ax.tick_params(axis='y',which='minor',left=False,labelleft=False)
    yticks=([df.exp0y.iloc[2]])
    yticklabels=[r'$\uparrow$'
                 '\n'
                 r'$y_P$'
                '\n'
                r'$\downarrow$']
    ax.set_yticks(yticks)
    ax.set_yticklabels(yticklabels)
    ax.plot(df.exp0x.iloc[1],df.exp0y.iloc[1],
           'o',color='black')
    ax.text(df.exp0x.iloc[1]-6,df.exp0y.iloc[1]+1,'I')
    ax.plot(df.exp0x.iloc[2],df.exp0y.iloc[2],
           'o',color='black')
    ax.text(df.exp0x.iloc[2]-6,df.exp0y.iloc[2]+20,'P')
    ax.axvline(df.exp0x.iloc[2],
               ymax=( math.log10(df.exp0y.iloc[2])-math.log10(ymin) ) \
               /    ( math.log10(ymax)            -math.log10(ymin) ),
               color='black',linestyle='--')
    ax.axhline(df.exp0y.iloc[2],
              xmax=df.exp0x.iloc[2]/xmax,
              color='black',linestyle='--')
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    ax.legend(loc='upper right')
    plt.text(0.03,0.85,'(c)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    plt.text(0.5,1.05,'Standard magnitude of the bloom peak',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    # --- EXP-3: SCHEMATIC
    ax=axs[4]
    h0=ax.plot(np.array(df.exp3x),np.array(df.exp3y),
               ':',color='violet',lw=3,label='EXP-3')
    ax.set_xlabel('Time')
    ax.xaxis.set_label_coords(0.05,-0.05)
    ax.set_xlim(xlims)
    xticks=([df.exp3x.iloc[2]])
    xticklabels=[r'$\leftarrow t_P \rightarrow$']
    ax.set_xticks(xticks)
    ax.set_xticklabels(xticklabels)
    ax.set_yscale('log')
    ax.set_ylabel('Log10(Chlorophyll $a$)')
    ax.yaxis.set_label_coords(-0.075,0.50)
    ymin=1E-4
    ymax=300
    ax.set_ylim(ymin,ymax)
    ax.tick_params(axis='y',which='minor',left=False,labelleft=False)
    yticks=([df.exp3y.iloc[2]])
    yticklabels=[r'$\uparrow$'
                 '\n'
                 r'$y_P$'
                '\n'
                r'$\downarrow$']
    ax.set_yticks(yticks)
    ax.set_yticklabels(yticklabels)
    ax.plot(df.exp3x.iloc[1],df.exp3y.iloc[1],
           'o',color='black')
    ax.text(df.exp3x.iloc[1]-6,df.exp3y.iloc[1]+1,
            'I\'\'\'',color='violet')
    ax.plot(df.exp3x.iloc[2],df.exp3y.iloc[2],
           'o',color='black')
    ax.text(df.exp3x.iloc[2]-6,df.exp3y.iloc[2]+20,
            'P\'\'\'',color='violet')
    ax.axvline(df.exp3x.iloc[2],
               ymax=( math.log10(df.exp3y.iloc[2])-math.log10(ymin) ) \
               /    ( math.log10(ymax)            -math.log10(ymin) ),
               color='black',linestyle='--')
    ax.axhline(df.exp3y.iloc[2],
              xmax=df.exp3x.iloc[2]/xmax,
              color='black',linestyle='--')
    rect=mpl.patches.Rectangle(xy=(xmin,ymin),
                               width=obs_breakup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
                               width=364-sim_freezeup,height=ymax,
                               linewidth=0,facecolor='grey',alpha=0.25)
    ax.add_patch(rect)
    ax.legend(loc='upper right')
    plt.text(0.03,0.85,'(d)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    plt.text(0.5,1.05,'Low magnitude of the bloom peak',
             horizontalalignment = 'center',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    # --- POSITION
    axs[0].set_position([0.13,0.92, 0.85,0.03])
    axs[1].set_position([0.13,0.66, 0.85,0.24])
    axs[2].set_position([0.13,0.64, 0.85,0.01])
    axs[3].set_position([0.13,0.34, 0.85,0.24])
    axs[4].set_position([0.13,0.04, 0.85,0.24])

with plt.style.context('mplstyles/exp.mplstyle'):    
    # Plot
    fig,axs=plt.subplots(5,1,figsize=(16,19.5))
    make_plots(axs)

    # --- SAVE
    plt.savefig('figures_progress/exp3.png',dpi=1000)


def make_plots(axs):
    xlims=(first_year[0],first_year[-1])
    xmin=xlims[0]
    xmax=xlims[1]
    
    # --- CONCEPTUAL MODEL
    ax=axs[0]
    im=PIL.Image.open('figures_progress/conceptual_model.png')
    ax.imshow(im)
    ax.get_xaxis().set_visible(False)
    ax.get_yaxis().set_visible(False)
    plt.text(0.03,0.86,'(a)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    # --- SCHEMATIC
    ax=axs[1]
    h0=ax.plot(np.array(df.exp0x),np.array(df.exp0y),
               '-',color='black',lw=3,label='EXP-0')
    ax.set_xlabel('Time')
    ax.xaxis.set_label_coords(0.05,-0.03)
    ax.set_xlim(xlims)
    xticks=([df.exp0x.iloc[1],df.exp0x.iloc[2]])
    xticklabels=[r'$\leftarrow t_I \rightarrow$',
                r'$\leftarrow t_P \rightarrow$']
    ax.set_xticks(xticks)
    ax.set_xticklabels(xticklabels)
    ax.set_yscale('log')
    ax.set_ylabel('Log10(Chlorophyll $a$)')
    ax.yaxis.set_label_coords(-0.075,0.50)
    ymin=1E-4
    ymax=300
    ax.set_ylim(ymin,ymax)
    ax.tick_params(axis='y',which='minor',left=False,labelleft=False)
    yticks=([df.exp0y.iloc[0],df.exp0y.iloc[2]])
    yticklabels=[r'$\uparrow$'
                 '\n'
                 r'$y_I$'
                '\n'
                r'$\downarrow$',
                r'$\uparrow$'
                 '\n'
                 r'$y_P$'
                '\n'
                r'$\downarrow$']
    ax.set_yticks(yticks)
    ax.set_yticklabels(yticklabels)
    ax.plot(df.exp0x.iloc[1],df.exp0y.iloc[1],
           'o',color='black')
    ax.text(df.exp0x.iloc[1]-6,df.exp0y.iloc[1]+1,'I')
    ax.plot(df.exp0x.iloc[2],df.exp0y.iloc[2],
           'o',color='black')
    ax.text(df.exp0x.iloc[2]-6,df.exp0y.iloc[2]+20,'P')
    ax.axvline(df.exp0x.iloc[1],
               ymax=( math.log10(df.exp0y.iloc[1])-math.log10(ymin) ) \
               /    ( math.log10(ymax)            -math.log10(ymin) ),
               color='black',linestyle='--')
    ax.axvline(df.exp0x.iloc[2],
               ymax=( math.log10(df.exp0y.iloc[2])-math.log10(ymin) ) \
               /    ( math.log10(ymax)            -math.log10(ymin) ),
               color='black',linestyle='--')
    ax.axvline(df.exp0x.iloc[2],
               ymax=( math.log10(df.exp0y.iloc[2])-math.log10(ymin) ) \
               /    ( math.log10(ymax)            -math.log10(ymin) ),
               color='black',linestyle='--')
    ax.axhline(df.exp0y.iloc[2],
              xmax=df.exp0x.iloc[2]/xmax,
              color='black',linestyle='--')
    ax.axhline(df.exp0y.iloc[0],
              xmin=df.exp0x.iloc[1]/xmax,
              xmax=df.exp0x.iloc[2]/xmax,
              color='black',linestyle='--')
    # ref.:
    # https://stackoverflow.com/questions/68111673/matplotlib-draw-ellipse-to-annotate-on-top-of-logarithmic-scale
    # answer by tmdavison
    # Arc centre coordinates
    x,y=df.exp0x.iloc[1],df.exp0y.iloc[1]
    # use the axis scale tform to figure out how far to translate
    arc_offset=mpl.transforms.ScaledTranslation(x,y,ax.transScale)
    # construct the composite tform
    arc_tform = arc_offset + ax.transLimits + ax.transAxes
    # Create the arc centred on the origin, apply the composite tform
    width,height=20,1 # not exact, by trial and error
    theta1,theta2=0,3 # not exact, by trial and error
    arc=mpl.patches.Arc(xy=(0,0),width=width,height=height,
                        color="grey",fill=False,lw=1,
                        theta1=theta1,theta2=theta2,
                        transform=arc_tform)
    ax.add_patch(arc)
    ax.text(174,1.2,'$m$')
    plt.text(0.03,0.86,'(b)',
             horizontalalignment = 'left',verticalalignment='center',
             transform=ax.transAxes, fontsize=25)
    
    # --- POSITION
    axs[0].set_position([0.13,0.50, 0.85,0.50])
    axs[1].set_position([0.13,0.04, 0.85,0.46])
    
with plt.style.context('mplstyles/exp.mplstyle'):    
    # Plot
    fig,axs=plt.subplots(2,1,figsize=(16,10.2))
    make_plots(axs)

    # --- SAVE
    plt.savefig('figures_progress/discussion.png',dpi=1000)

	depth	chlHPLC	doy
20357	0.5	0.1056	146
20358	1.5	0.0899	146
20359	5.0	0.0906	146
20360	10.0	0.0905	146
20361	20.0	0.1011	146
20362	40.0	0.0633	146
20363	0.5	0.0971	148
20364	1.5	0.1483	146
20365	5.0	0.1086	148
20366	10.0	0.0890	146
20367	20.0	0.0735	148
20368	40.0	0.0399	148

	depth	chlHPLC	doy
20357	0.5	0.1056	146
20358	1.5	0.0899	146
20359	5.0	0.0906	146
20360	10.0	0.0905	146
20361	20.0	0.1011	146
20362	40.0	0.0633	146
20363	0.5	0.0971	148
20364	1.5	0.1483	148
20365	5.0	0.1086	148
20366	10.0	0.0890	148
20367	20.0	0.0735	148
20368	40.0	0.0399	148

Introduction¶

Directory structure¶

Environment¶

First, load libraries¶

Maps¶

Maps: stations¶

Maps: plot¶

Select year¶

Grid¶

Observations¶

Observations: ice concentration¶

Observations: snow and ice thickness¶

Observations: snow and ice thickness: load dataset¶

Observations: snow and ice thickness: select¶

Observations: PAR¶

PAR¶

Observations: PAR: Load dataset¶

Observations: PAR: Select¶

Observations: PAR: calculate¶

Observations: PAR: 3 m¶

Observations: PAR: isolume¶

Observations: nutrients¶

Observations: nutrients: load dataset¶

Observations: nutrients: select¶

Observations: nutrients: 3 m¶

Observations: total dissolved inorganic iron¶

Observations: total dissolved inorganic iron: load dataset¶

Observations: chlorophyll a by HPLC¶

Observations: chlorophyll a by HPLC: load dataset¶

Observations: chlorophyll a by HPLC: select¶

Observations: chlorophyll a by HPLC: vertically integrated chlorophyll a (0-40m)¶

Observations: chlorophyll a by HPLC: vertically integrated chlorophyll a (0-100m)¶

Observations from literature¶

Observations from literature: maximum photosynthesis rate at saturation irradiance in nutrient replete conditions ($P^C_{max} \gamma^T$, $\mathrm{ d^{-1} }$)¶

Observations from literature: lineare initial slope of the Chl a-specific photosynthesis versus irradiance curve ($\alpha^{chl}$, $\mathrm{ g\ C\ (g\ Chl\ a)^{-1}\ h^{-1}\ (\mu mol\ photons\ m^{-2}\ s^{-1})^{-1} }$)¶

Observations from literature: light saturation parameter ($E_k$, $\mathrm{ \mu mol\ photons\ m^{-2}\ s^{-1} }$)¶

Outputs from model NEMO-LIM3¶

Output from model NEMO-LIM3 and input into package gud of MITgcm: sea ice concentration¶

Output from model NEMO-LIM3 : snow volume¶

Output from model NEMO-LIM3 : ice volume¶

Model¶

Model: traits¶

Model: traits: set indir¶

Model: traits: create RunInfo object¶

Model: traits: Equivalent spherical diameters (ESD, $\mathrm{ \mu m }$)¶

Model: traits: Maximum growth rate ($P^C_{max}$, $\mathrm{ d^{-1} }$)¶

Half saturation for growth on nitrate ($k_{NO_3}$, $\mathrm{ mmol\ N\ m^{-3} }$)¶

Maximum grazing rate ($g_{max}$, $\mathrm{ d^{-1} }$)¶

Model: outputs¶

Model: PAR¶

Model: mean PAR: 3 m¶

Model: mean PAR: isolume $\mathrm{ z_{0.415}=0.415\ mol\ photons\ m^{-2} d^{-1} }$¶

Model: nutrients¶

Model: nutrients: 3 m¶

Model: total dissolved inorganic iron¶

Model: chlorophyll a¶

Model: chlorophyll a: vertically integrated chlorophyll a (0-40m)¶

Model: chlorophyll a: vertically integrated chlorophyll a (0-100m)¶

Model: chlorophyll $a$ in simulation EXP-1¶

Model: chlorophyll a in simulation EXP-1: vertically integrated chlorophyll a (0-100m)¶

Model: chlorophyll $a$ in simulation EXP-2¶

Model: chlorophyll a in simulation EXP-2: vertically integrated chlorophyll a (0-100m)¶

Model: chlorophyll $a$ in simulation EXP-3¶

Model: chlorophyll a in simulation EXP-3: vertically integrated chlorophyll a (0-100m)¶

Model: temperature¶

Model: $\gamma^T$¶

Model: maximum photosynthesis rate at saturation irradiance in nutrient replete conditions ($P^C_{max} \gamma^T$, $\mathrm{ d^{-1} }$)¶

Model: lineare initial slope of the Chl a-specific photosynthesis versus irradiance curve ($\alpha^{chl}$, $\mathrm{ g\ C\ (g\ Chl\ a)^{-1}\ h^{-1}\ (\mu mol\ photons\ m^{-2}\ s^{-1})^{-1} }$)¶

Model: light saturation parameter when there is light or no light ($E_k$, $\mathrm{ \mu mol\ photons\ m^{-2}\ s^{-1} }$)¶

Model: light saturation parameter when there is light ($E_k$, $\mathrm{ \mu mol\ photons\ m^{-2}\ s^{-1} }$)¶

Plot¶

Plot: sea ice, PAR, nutrients, Chl a, and $z_{0.415}$¶

Plot: traits¶

Plot: nitrate, silicic acid, phosphate, iron, and chlorophyll $a$¶

Plot: compare traits from literature data and simulation¶

Plot: EXP-1 (no light under sea ice)¶

Plot: EXP-2 (no minimum biomass)¶

Plot: EXP-3 (differing wintertime nitrate)¶

Plot: discussion¶